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general  viscous  damping  model  is  used  to  verify  the  calculated  modal  parameters. 
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I. 


INTRODUCTION 


Mechanical  energy  transforms  into  heat  and  dissipates  in  all  vibrating  systems. 
There  are  many  energy  dissipation  mechanisms  that  contribute  to  the  damping  in  the 
structure  system.  Some  of  these  mechanisms  are:  fluid  resistance  and  coupling,  internal 
friction  (material  damping),  and  friction  at  a  joint.  All  of  these  dissipation  mechanisms 
have  been  shown  to  be  a  function  of  many  variables,  including  a  structure’s  shape  or 
geometry,  its  material  properties,  temperature,  frequency,  boundary  conditions,  and 
different  excitation  energy  levels.  Usually  over  90  percent  of  the  inherent  damping 
associated  with  fabricated  build-up  structures  originates  in  the  mechanical  joints  (Beards 
and  Woodwat,  1985).  These  mechanical  joints  are  friction  joints,  which  dissipate  energy 
during  the  vibration  of  a  structure.  Reducing  the  contact  force  in  bolted  structural 
connections  can  reduce  system  vibration  amplitudes  by  enhancing  joint  damping  capacity 
(Shin  et  al.,  1991). 


Naval  ship  structure  systems  have  mostly  welded  joints  and  all  stiffeners  are  also 
welded  to  hull  plates,  decks  and  bulkheads.  The  ship  system  also  has  many  energy 
dissipation  sources  such  as  long  cable  trays,  hangers,  snubbers,  the  surrounding  fluid 
coupled  with  ship  hull,  etc.  The  ship  system  damping  is  measurable,  but  difficult  to 
quantify  (Rutgerson,  2002).  In  conjunction  with  ship-shock  simulation  based  transient 
analysis,  time-domain  representation  of  system  damping  is  desirable  using  the  frequency- 
domain  characteristics  of  damping.  The  damping  studies  used  for  analysis  were 
conducted  using  2  sec  data  from  the  DDG  53  Ship  Shock  Trials. 


The  goal  of  this  study  is  to  present  the  damping  model  in  Rayleigh  damping  form 
of  a  naval  ship  system  for  ship  shock  transient  time  domain  analysis.  In  this  study,  the 
Complex  Exponential  Method  is  used  for  extraction  of  modal  parameters  in  the  time 
domain.  The  Inverse  Fourier  Transform  of  Mobility  form  of  the  general  viscous  damping 
model  verifies  the  calculated  modal  parameters.  Two  factors  in  the  Rayleigh  damping 


1 


model  are  calculated  using  modal  frequency  and  modal  damping  ratios.  The  statistical 
characteristics  of  two  Rayleigh  factors  are  quantified  in  each  categorized  area.  Then  the 
spatially  dependent  Rayleigh  damping  model  is  investigated  and  a  model  to  be  used  in 
shock  transient  analysis  is  recommended. 
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II.  THEORY  OF  GENERAL  VISCOUSLY  DAMPED  SYSTEM 


The  general  equation  of  motion  for  a  MDOF  (Multi  Degree  Of  Freedom)  system 
with  viscous  damping  and  harmonic  excitation  is: 


[M]{x}  +  [C]{i}  +  [X]{x}={/}  (l) 


In  the  above  equation,  [M]  is  the  system  mass  matrix,  [K]  is  the  system  stiffness 
matrix,  [C]  is  the  system  damping  matrix,  {x}  is  the  system  response  vector  and  {f}  is  the 
forcing  vector. 

We  consider  first  the  case  where  there  is  zero  excitation  in  order  to  determine  the 
natural  modes  of  the  system  and  to  this  end,  we  assume  a  solution  to  the  equations  of 
motion  which  has  the  form 


(2) 


Substituting  this  into  the  appropriate  equation  of  motion  gives: 

(MM  +  S  [c]+[*])  M  =  {0}  (3) 


the  solution  of  which  constitutes  a  complex  eigenvalue  problem.  In  this  case, 
there  are  2N  eigen  values,  sr  in  complex  conjugate  pairs.  (This  is  an  inevitable  result  of 

the  fact  that  all  the  coefficients  in  the  matrices  are  real  and  thus  any  characteristic  values, 
or  roots,  must  either  be  real  or  occur  in  complex  conjugate  pairs.)  There  is  an 
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eigenvector  corresponding  to  each  of  these  eigenvalues,  but  these  also  occur  as  complex 
conjugates.  Hence  we  can  describe  the  eigensolution  as: 


Sr’  S*  and  ,  {S*}r  r  =  1,  N 


(4) 


It  is  customary  to  express  each  eigenvalues  sr  in  the  form 


s=a>r(-^r+i^C)  (5) 


where  cor  is  the  ‘natural  frequency’  and  is  the  critical  damping  ratio  for  that 
mode.  Sometimes,  the  quantity  cor  is  referred  to  as  the  ‘undamped  natural  frequency’  but 

this  is  not  strictly  correct,  except  in  the  case  of  proportional  damping  (or,  of  course,  of  a 
single  degree  of  freedom  system). 

The  eigen  solution  possesses  orthogonality.  In  order  to  examine  these  we  must 
first  note  that  any  eigenvalue/eigenvector  pair  satisfies  the  equation 


(5)[M]+.,.[C]+[X])V}r={0}  (6) 


and  then  we  pre-multiply  this  equation  by  {y/}Tq  so  that  we  have: 


Wl(^M+^[c]+[x])Mr  =  {o}  (7) 
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A  similar  expression  to  (6)  can  be  produced  by  using  A  and  {y/}q : 


«M+^[c]+[x]){r}9={o}  (8) 

which  can  be  transposed,  taking  account  of  the  symmetry  of  the  system  matrices, 

to  give: 

Ml(<[MK[c]+[x])  =  {o}r  (9> 

If  we  now  post-multiply  this  expression  by  {y/}r  and  subtract  the  result  from  that 
in  Equation  (7),  we  obtain: 

h =o  do) 

and  provided  sr  and  sq  are  different,  this  leads  to  the  first  of  a  pair  of 
orthogonality  equations: 

(*-•+*,  )Ml  MM, + M^MM,  =  °  ("> 

A  second  equation  can  be  derived  from  the  above  expressions  as  follows: 

Multiply  (7)  by  s  and  (9)  by  sr  and  subtract  one  from  the  other  to  obtain: 
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(12) 


v,  Ml  MM,  -  Ml  [K\W)r  =  o 


These  two  equations  -  (11)  and  (12)  -  constitute  the  orthogonality  conditions  of 
the  system  and  it  is  immediately  clear  that  they  are  far  less  simple.  However,  it  is 
interesting  to  examine  the  form  they  take  when  the  modes  r  and  q  are  found  as  a  complex 
conjugate  pair.  In  this  case,  we  have  that 


S,  =©,(-£.  -/-y/  1-M  (13) 


and  also  that, 


=M), 


(14) 


Inserting  these  into  Equation  (11)  gives 


-2  (0r£r  M}(  [AH  {v},  +  M}(  [C]  {H, =0  (15) 


from  which  we  obtain: 


MlMiM  c, 

{v'frWW},  ™r 


(16) 


Similarly,  inserting  (13)  and  (14)  into  (12)  gives 
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^{Y}TrWW}r  -  {Y}TAKW}r  =  0 


(17) 


from  which 


{yJXBMz  k 

{v'frWW},  m 


In  these  expressions,  mr ,  kr ,  and  cr  may  be  described  as  modal  mass,  stiffness 

and  damping  parameters  respectively  although  the  meaning  is  slightly  different  to  that 
used  in  the  other  systems. 


A.  FORCED  RESPONSE  ANALYSIS 


Returning  to  Equation  (1),  and  assuming  a  harmonic  response: 


W01 


(19) 


we  can  write  the  forced  response  solution  directly  as 


{x}  =[[X]  -  ffl2[M]  +  iffl[C]  ]"'  {/}  (20) 
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but  this  expression  is  not  particularly  convenient  for  numerical  application. 
Define  a  new  coordinate  vector  {y},  which  is  of  order  2N,  and  which  contains  both  the 
displacements  {x}  and  the  velocities  {ij  : 


{y}  = 


(21) 


Equation  (1)  can  then  be  written  as: 


[C :  M]Nx2N  {j>}2Kx1  +  [K :  0]  {y}  =  {0}M  (22) 


However,  in  this  form  we  have  N  equations  and  2N  unknowns  and  so  we  add  an 
identity  equation  of  the  type: 


[M  :  0]  {y}  +  [0 :  -M]  {y}  =  {0}  (23) 


which  can  be  combined  to  form  a  set  of  2N  equations 


fC 

KM 


{>}  + 


(K 

v0 


0  A 

-Mj 


=  (0} 


(24) 


which  can  be  simplified  to: 
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[-4]  W  +  [S]  W  =  {0} 


(25) 


These  equations  are  now  in  a  standard  eigenvalue  form  and  by  assuming  a  trial 
solution  of  the  form  {y}  =  {y}es> ,  we  can  obtain  the  2N  eigenvalues  and  eigenvectors  of 

the  system,  2.  and  {$},. ,  which  together  satisfy  the  general  equation: 


(\[A]  +  [B]){6}r  =  { 0};  r  =  1, 2N  (26) 


These  eigen  properties  will,  in  general,  be  complex  although  for  the  same  reasons 
as  previously  they  will  always  occur  in  conjugate  pairs.  They  possess  orthogonality 
properties,  which  are  simply  stated  as 


{ovum 


• a 


r. 


{oYiBm 


(27) 


and  which  have  the  usual  characteristic  that 


r  =  1 ,  2N 


(28) 


Now  we  may  express  the  forcing  vector  in  terms  of  the  new  coordinate  system  as: 
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(29) 


{P} 


2Nxl 


and  assuming  a  similarly  harmonic  response  and  making  use  of  the  previous 
development  of  a  series  form  expression  of  the  response.  We  may  write: 


r= i  ar(ico-sr) 


(30) 


However,  because  the  eigenvalues  and  eigenvectors  occur  in  complex  conjugate 
pairs,  this  last  equation  may  be  written  as: 


N  f{0}TAp}{e}r  ,  {£*}>} 


=  Z 

r= 1 


+ 


ar  (ico  —  s)  ar  (ico  —  s) 


(31) 


At  this  stage,  it  is  convenient  to  extract  a  single  response  parameter,  say  , 
resulting  from  a  single  force  such  as  fk  -  the  receptance  frequency  response  function,  a1’ , 
and  in  this  case  Equation  (31)  leads  to: 


X,. 


N 


( 


-L  =  a)(m)=Yd 

Jk  r= 1 


6.  6. 

r  j  r  k 


+ 


e  e: 

r  j  r  k 


Kar(ico-sr )  ar{ico-sr)J 


(32) 


or, 
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N 


a 


»=E 


(  A  A*  ^ 

r  Ajk  r  Ajk 


r- 1 


(ia)-sr)  ( ico-s*r)j 


(34) 


where,  ^;/c  = 


r9j  A 

a,. 


A  - 

r^jk  ~ 


n  a 

r  “j  ftk 

* 

a,. 


Equation  (34)  describes  the  displacement  response  at  ‘j’  degree  of  freedom 
under  excitation  at  ‘k’  degree  of  freedom  of  general  viscously  damped  system  in 
frequency  domain.  If  Fourier  Transform  is  made  onto  Equation  (34),  we  can  get  the 
Impulse  Response  Function  hjk  (t)  of  this  system. 


i  00  N 


f 


t-Ajk 


A 

r^jk 


\ 


V 


+ 

( i(D-sr )  (ico  -  s'. )  j 


icot 


dco 


(35) 


If  we  let  ( ico-sr )  as  iz ,  then,  by  Cauchy’s  Integral  Formula,  it  becomes: 


-|  uu 

2- 1 

9  7 r  J 


00  A 

r  jk  Jcot 


2  7ri(io-sr) 


e  da>  = 


A  s/A  oo  jzt 
r  J* 


2m  v  z 

— oo 


J  — dz  =  rA 


-jk^  (36) 


and  similarly, 


1  00  A 

J_  f  '  AJk 

0  7r  J  ( irn  —  < 


2/r  ^  (ico  -  sr )  2 m 


A‘ev  "  -» 


— }  —dz  =  rA’lles:'  (37) 
z 
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From  Equations  (36)  and  (37),  the  Impulse  Response  Function  of  a  general 
viscous  damping  system  can  be  written  as: 


h 


jk 


N  / 

<o=Z( 

r= 1 


A 


sJ 


~jk 


+  A 


sJ 


(38) 


The  forced  response  can  then  be  calculated  as  Equation  (39)  in  time  domain. 


t 

Xj=\fk(T)hjk(t-T)dT  (39) 

0 
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III.  PROCEDURE  OF  DAMPING  CALCULATION  FROM 

MEASURED  DATA 

A.  MODAL  PARAMETER  EXTRACTION 

In  Equation  (34),  r  Ajk ,  rA*jk,  sr,  s*r  are  called  modal  parameters, 

r  Ajk ,  rA*  are  called  modal  constants,  and  eigenvalues  sr,  s*r  contain  information  of 

modal  properties  such  as  modal  frequency  cor  and  modal  damping  ratio  Cr .  Calculation 

of  damping  ratio  from  the  measured  shock  trial  data  needs  utilizing  modal  parameter 
extraction  methods,  the  simple  3-dB  (half  power)  bandwidth  measurement  or  logarithmic 
decay  rate  calculation  can  be  incorrect  in  real  cases,  because,  in  usual  cases,  real 
measurement  data  not  only  contain  noise  components  but  also  have  many  closely  coupled 
frequency  components. 

Basically  there  are  two  groups  of  techniques  in  the  field  of  experimental  modal 
analysis.  One  is  related  to  frequency  domain  analysis  methods  that  use  Frequency 
Response  Function  of  measured  input  and  output  data.  This  group  of  methods  is  widely 
used,  from  single  degree  of  freedom  circle  fitting  to  complex  multi-degree  of  freedom 
fitting  methods.  And  the  other  is  related  to  time  domain  analysis  techniques  that  use 
Impulse  Response  Function  as  analysis  data. 

In  this  study,  the  time  domain  method  is  used  because  the  measured  useful 
shock  trial  data  sets  are  too  short  in  time  to  obtain  sufficient  frequency  resolution  in  the 
frequency  domain,  and  there  is  only  one  event,  shot  data,  thus  the  averaging  process 
cannot  be  done.  Complex  Exponential  Method  (CEM),  one  of  the  effective  modal 
parameter  extraction  methods  in  the  time  domain,  is  used  to  extract  modal  parameters. 

The  2  second  long  measured  data  has  been  used  in  this  study  and  the  effective 
frequency  span  is  limited  from  3  Hz  to  250  Hz.  The  original  measured  data  has  been 
band-pass  filtered  from  2  Hz  to  250  Hz  to  avoid  long-term  trends  in  low  frequencies  and 
to  remove  unwanted  high  frequency  noise  components. 


14 


1.  Complex  Exponential  Method  (CEM) 


In  the  field  of  experimental  modal  analysis,  a  term  Receptance  is  widely  used  to 
describe  the  ratio  of  displacement  response  to  excitation  force,  and  the  term  Mobility  is 
used  in  describing  the  ratio  of  velocity  response  to  excitation  force.  The  Receptance 
(displacement/force)  a(a>)  of  a  general  viscously  damped  system,  Equation  (34),  can  be 

rewritten  as  Equation  (40),  complex  eigenvalue  sr  is  as  in  Equation  (40. a), 


.=1  J<o-sr  jco-s 


(40.  a) 


or, 


2N  A 

a(oj)  =  V -  ;  sr^s*r,  4^4*  ,for  r>N  (40.b) 

r=xjco-sr 


we  can  get  velocity  v  =  Vej<M  by  time  differentiating  displacement  x  =  XejM ,  that 
is, 


v=V(co)eim  =jmXejm  () 


and  Mobility( velocity/force)  Y(oj)  can  be  related  to  a(co)  (displacement/force), 
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Y(co)  =  jcoa{co) 


(41) 


The  corresponding  Impulse  Response  Function  (IRF)  can  be  obtained  by 
taking  Inverse  Fourier  Transform  (IFT)  of  the  Receptance  a(co)  as  Equations  (35) 
through  (38), 


IN 

h(t)  =  YJ  A^'  C42) 

r- 1 

By  time  differentiating,  the  velocity  form  of  IRF  can  be  expressed  as, 

2N 

h  (t)  =  Yj  A^N'1  («) 

r= 1 


Then  the  measured  velocity  data  set  sampled  by  D  t  time  steps  can  be  expressed 
as  follows: 


/?,,  A2, . ,  hq  =h( 0),  h(At),  h(2At), . ,  h(qAt ) 


(44) 


Hereafter,  the  exponential  term  will  be  simplified  using  the  following  notation. 


esAt 


->K 


(45) 


Thus  for  the  j-th  sample  data,  the  Equation  (43)  becomes, 
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(46) 


2N 


hi=Y,ArSrVj 

r— 1 


which,  when  extended  to  the  full  data  set  of  q  samples  (j=l,2,....,q),  gives: 


II 

-sp 

6)  A 

+ 

S2A2  +  •••• 

'  *  +  S2nAn 

K  = 

A 

+ 

V2  s2A2  h — 

* "  +  AnS2nAn 

II 

hV, 

+ 

V2  s2A2  h — 

+  V2N  s2NA2N 

II 

•-sP 

+ 

V2qs2A2  +•••• 

"  ‘  +  AA  S2nAn 

Provided  that  the  number  of  sample  points  q  exceeds  4N,  this  equation  can  be 
used  to  set  up  an  eigenvalue  problem,  the  solution  yields  the  complex  natural  frequencies 
contained  in  the  parameters  VX,V2,  etc. 

Multiply  each  equation  in  (47)  by  a  coefficient,  J3,  to  form  the  following  set  of 
equations: 


II 

-s? 

< 

PoA  + 

PoA 2  + 

•  ••  +  P§a2n 

PA  — 

pm + 

PM2  +■••• 

•  •  +  P\AnAn 

11 

OS  ” 

M24  + 

P2v22a2  +•• 

+  PA  in  An  (48) 

II 

•  ^  ^ 

Pqv?\  + 

P/2"a2+- 

+  PqAA  An 

Adding  all  equations  in  (48)  vertically  results  in, 
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(49) 


q  2N  q 

zM=zmw) 

i=0  7=1  i=0 

The  coefficients  [3J  s  are  taken  to  be  the  coefficients  in  the  polynomial  equation, 

/?0+/?1r+/?2r2+/?3r3  + . +/?,r=o  m 

The  roots  are  Vt,  V2,  Vq  . 

Next,  the  values  of  the  (1  coefficients  are  to  be  sought  in  order  to  determine  the 
roots  of  Equation  (50)  -  values  of  Vr  -  and  hence  the  system  natural  frequencies.  Now, 

recalling  that  q  is  the  number  of  degrees  of  freedom  of  the  system  model.  It  is  now 
convenient  to  set  these  two  parameters  to  the  same  value,  i.e.  letq  =  2N . 

If  we  find  P  coefficients  that  make  Equation  (50)  fulfilled,  then  Equation  (50) 
can  be  expressed  as, 

IN 

Y^PjVj  =0  ;  r  =  1,  2N  (51) 

7=0 

And  thus  every  term  on  the  right-hand  side  of  Equation  (49)  is  zero. 
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(52) 


IN 


XM=  o 


Rearranging  Equation  (52),  by  moving  the  last  term  of  left-hand  side  to  right-hand 

side, 


2N-\ 

^  A  4  =  —  \n 


1=0 


by  setting  =  1 


(53) 


Repeat  the  process  from  (44)  to  (53)  using  different  set  of  IRF  data  points  and 
further  choose  the  new  data  set  that  overlaps  considerably  with  the  first  set  -  In  fact,  for 
all  but  one  item. 


Successive  applications  of  this  procedure  lead  to  a  full  set  of  2N  equations: 


K 

K 

h2 

^2N-l 

'  A  ' 

C  * 

h 2N 

K 

ii  2 

•••  j**- 

4ziV 

< 

A 

>=  —  < 

h 

"2jV+1 

> 

1 - 

fo 

\n 

^2N+\ 

1 

<N 

1 

^2N-\ , 

pAN-X, 

-^-2Nx2N  ^2Nxl  W  2Nx1 
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The  unknown  coefficients  {0}  can  be  found  from  Equation  (54).  Now  the 
values  of  V] ,  V2,  ...,  V2N  can  be  determined  using  Equation  (50)  and  subsequently  the 
system  natural  frequencies  can  be  found  using  the  following  relationship. 


V=e 


s„At 


(55) 


Using  Equation  (47),  corresponding  modal  constants^,  A2 ,  ...,  A1N  can  be 
calculated,  this  may  be  written  as, 


'  1 

1 

i 

1 

As  i 

K 

h 

v2 

y3 

v 

v2  N 

A^2 

K 

V2 

f2 

v2 

■  V  2 

v2  N 

... 

< 

AS3 

>  =  < 

h2 

1 

^  . 

fo 

v2N~' 

V2N-, 

fo 

1 _ 

AlNS2N „ 

jh  n-u 

or 

[V]  {A}  =  {h}  (57) 
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2. 


Verification  of  Extracted  Modal  Parameters 


The  modal  parameters  calculated  according  to  the  mentioned  procedure  can  be 
verified  by  comparing  synthesized  time  histories  to  the  originally  measured  time 
histories.  Mobilities  can  be  determined  by  frequency  differentiation,  multiplied  by  jco , 
from  Equations  (40)  and  (41), 


ha>)  =  £  ;  sr  =  -corfr  +  jcor 

r=lJCO-Sr  JCO-Sr 


(58) 


or 


f(<y)  =  £-^-  ;  Sr  =>  s',  Ar  =>  4  ,  for  r>  N  (59) 

r=lJCO-Sr 

In  Equations  (58)  and  (59),  extracted  modal  parameters  Ar ,  A* ,  sr ,  s*r  are  used 

to  calculate  the  synthesized  Frequency  Response  Function  in  frequency  domain.  And  by 
inverse  Fast  Fourier  Transform  (IFFT),  taking  the  real  part  of  the  results,  synthesized  IRF 
can  be  calculated  and  compared  to  original  time  histories.  Frequency  bandwidth  A/ 
should  be  multiplied  during  calculation  procedure  to  generate  band  level. 


M  =  Real  Part  of  Inverse  FFT ( At  f ) 


(60) 
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3. 


Calculation  of  Rayleigh  Damping 


Rayleigh  damping  is  a  kind  of  general  proportional  damping  model.  It  assumes 
that  damping  matrix  [C]  in  Equation  (1)  can  be  represented  as  linear  combination  of  the 
mass  matrix  and  stiffness  matrix.  Then  the  damping  matrix  can  be  easily  decoupled  to  the 
modal  damping  matrix.  Using  the  Rayleigh  damping  representation,  the  damping  matrix 
can  be  represented  as, 


[C]=a[M]+p[K 


(61) 


or  by  using  mass  normalized  modal  matrix  [<^>] , 


>f  [C]M=[2®^]^=«/ +P 


(62) 


By  using  Equation  (62)  for  all  2N  modes,  following  2N  equations  can  be  set. 


a+pa\  -2  (d£x 

a+  Bcor,  =2 

1  z  z  (63) 

a + Pcq2N  =  2co2N^2n 
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If  2N  is  larger  than  2,  then  Equation  (64)  becomes  over-determined,  with  2 
unknowns  and  2N  equations.  By  post-multiplying  the  transpose  matrix  of  [W]  to  both 
sides  of  Equation  (64),  we  can  get  Equation  (65). 


2x2  N 


(65) 


Then  two  Rayleigh  parameters  a  and  (1  are  calculated  as, 


=  (\W]T  ^[W\NJ\W]T  2«2„ 


2iVxl 


(66) 


As  a  final  step,  modal  damping  ratios  for  each  vibration  modes  are  calculated  in 
each  categorized  area  of  ship. 


1 

a  — 

2  co 


(67) 
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IV.  RESULTS  OF  MODAL  PARAMETER  EXTRACTION 


A.  VERIFICATION  RESULTS 


Modal  parameters  for  a  total  of  773  sets  of  measured  data  are  calculated  for  each 
measuring  position  and  direction. 

The  following  figures  show  some  results  of  the  modal  parameter  extraction 
method  according  to  the  aforementioned  procedure.  Figure  1,  Figure  3,  and  Figure  5 
show  the  original  measured  data  sets,  while  Figure  2,  Figure  4,  and  Figure  6  show 
synthesized  curves  with  the  calculated  modal  parameters  cor ,  t,r  and  Ar .  The  black  line 

shows  parts  of  the  original  signal  between  the  125  msec  point  and  the  1200  msec  point. 
The  red  line  stands  for  the  synthesized  curve  of  those  parts.  The  125  msec  parts  of  each 
of  the  measured  data  sets  are  not  included  in  analysis  to  avoid  being  mixed  with  the  effect 
of  excitation  signal.  Also  the  latter  parts  that  contain  the  secondary  excitation  are  not 
included  in  the  analysis.  This  secondary  excitation  can  be  seen  on  Figure  1,  Figure  3,  and 
Figure  5  around  the  1250  msec  point. 


F igure  1 .  Measured  Data  at  A200 1 AI 
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Figure  2.  Synthesized  Results  at  A2001A  (between  125  msec  and  1200  msec) 


Figure  3.  Measured  Data  at  A2004A 


Figure  4.  Synthesized  Results  at  A2004A  (between  125  msec  to  1200msec) 
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Figure  5.  Measured  Data  at  A3  5  06V 


Time  (m&ee) 


Figure  6.  Synthesized  Results  at  A3506  (between  125  msec  and  1200  msec) 


B.  RESULTS  OF  RAYLEIGH  DAMPING  CALCULATION 


A  total  of  773  data  sets  were  categorized  into  67  area  groups,  based  on  the 
location  of  the  measuring  sensor  installation,  as  well  as  by  the  direction  of  measurement, 
which  were  the  athwartship  and  vertical  directions  of  the  ship. 

Figure  7  through  Figure  10  show  some  typical  results  of  the  Rayleigh 
Damping  coefficients  calculated  according  to  the  aforementioned  procedure,  using 
Equations  (61)  through  (67).  Each  black  square  point  is  a  mode.  The  red  lines  represent 
the  regenerated  modal  damping  ratio  using  a  and  (3  calculated  in  least-square  sense. 
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Each  figure  tagged  as  ‘Original’  in  caption  is  from  the  initial  step  results,  whereas  the 
others  are  of  the  final  results,  which  have  been  modified  by  eliminating  the  unrealistic 
and  noisy  components.  From  the  ‘Original’  figures,  we  can  identify  that  some  of  the 
results  are  scattered  and  contain  damping  ratios,  which  cannot  be  regarded  as  reasonable 
against  the  physical  sense.  Some  modifications  to  the  above  results  have  been  made.  The 
unreasonable  damping  data  points  have  been  removed  from  curve  fitting  a  and  p.  The 
modes  of  which  modal  constants  Ar  are  seriously  less  than  one  thousandth  of  the 

maximum  value  in  each  measuring  position  have  been  removed.  The  modes  that  contain 
damping  ratios,  which  are  greater  than  0.5  have  been  removed.  Likewise,  the  modes  with 
damping  ratios  that  contain  a  great  deal  of  scattered  from  the  initial  curve-fitted  have 
been  also  removed.  The  final  results  are  shown  in  both  of  linear  scales  and  logarithmic 
scales. 


The  figures  presented  in  Appendix  A  illustrate  the  damping  calculation  results  for 
the  remaining  areas  that  were  studied. 


Figure  7.  Modal  Damping  Ratio  at  Area  6,  Athwartship  Direction  (Original) 
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Damping  Ratio  Damping  Ratio 


(in  Linear  Scale) 


(in  Logarithmic  Scale) 


Figure  8.  Modal  Damping  Ratio  at  Area  6,  Athwartship  Direction  (Modified) 


Figure  9. 
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Figure  10.  Modal  Damping  Ratio  at  Area  6,  Vertical  Direction  (Modified) 
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C.  CURVE-FITTED  RAYLEIGH  DAMPING  a  AND  p  FOR  EACH  AREA 


The  Curve-fitted  Rayleigh  damping  coefficients,  a  and  p,  for  each  categorized 
area  are  presented  in  Table  1  and  Table  2,  for  the  athwartship  and  vertical  directions, 
respectively.  Figure  2  is  a  profile  view  of  the  DDG  51  Arleigh  Burke  Class  Destroyer. 
This  drawing  shows  the  major  transverse  frame  positions  of  the  ship. 


Figure  1 1 .  Transverse  Frame  Locations  of  the  DDG  5 1  Class  Destroyer 


T able  1 .  Rayleigh  Damping  Results  for  the  Athwartship  Direction  by  Area  Group 


Athwartship 
Position 
p  =  port 
s  =  starboard 


Tp-Ts 


CL-5's 


20’p-20’s 


17's 


CL-15's 


5’p-CL 


_ 8'P _ 


5’p-5’s 


25’p-20’p 


lTp-ll’s 


3Ep-3Es 


22’p-CL 


6's-lO’s 


lO’p-ll’s 


2.09E+01 


1.91E+01 


Alpha 


2.05E+01  1.22E-06 


2.14E+01  4.69E-06 


1.91E+01  5.08E-07 


1.29E+01  1.11E-05 


L.90E-06 


1.79E+01  1.04E-06 


1.88E+01  5.74E-06 


1.71E+01  4.23E-06 


1.90E+01  3.07E-06 


1.76E+01  1.34E-06 


2.04E+01  2.44E-07 


L.72E-06 


Number  of 
modes 
used  in 
curvefitting 
a  and  (3 


2.03E+01  2.1  IE-06 


1.95E+01  2.49E-07 
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Area 

No. 

Deck 

Frame 

Athwartship 
Position 
p  =  port 
s  =  starboard 

Alpha 

Beta 

Number  of 
modes 
used  in 
curvefitting 
a  and  P 

31 

2.5 

398 

CL-8's 

1.67E+01 

1.50E-05 

11 

32 

3 

192.5 

6'p-23's 

2.14E+01 

2.63E-06 

26 

33 

3 

277 

20’p-14’s 

1.98E+01 

4.04E-06 

53 

35 

3 

129 

6’p-CL 

2.45E+01 

5.90E-06 

11 

37 

3 

151 

36’p-36’s 

1.76E+01 

2.34E-06 

59 

40 

3 

331.5 

3’s-CL 

1.84E+01 

2.18E-05 

11 

41 

3 

397.5 

13'p-CL 

1.94E+01 

1.39E-06 

54 

43 

3.5 

156 

8'p-2’s 

2.00E+01 

8.14E-07 

125 

48 

4 

174 

1  l'p-1  l’s 

1.35E+01 

1.91E-06 

21 

49 

4 

232.5 

6'p-CL 

2.23E+01 

-1.26E-08 

33 

53 

5 

280 

CL 

1.66E+01 

1.16E-05 

19 

54 

5 

322.5 

17'p-15's 

1.90E+01 

1.30E-06 

75 

55 

6 

150 

CL 

2.05E+01 

3.75E-06 

29 

56 

6 

160 

CL 

1.97E+01 

3.11E-06 

21 

57 

7 

150 

CL 

1.77E+01 

1.30E-05 

8 

58 

7.5 

155.5 

8's-CL 

1.11E+01 

8.42E-06 

22 

60 

8 

173.5 

32'p-32's 

1.40E+01 

4.19E-06 

63 

61 

9 

168 

CL 

1.17E+01 

2.69E-06 

32 

62 

9 

171.5 

25'p-25’s 

1.27E+01 

3.82E-06 

33 

64 

10 

174 

6’s-lO’s 

1.93E+01 

3.48E-06 

17 

65 

11 

181 

CL 

1.58E+01 

7.18E-06 

12 

67 

MD 

281.5 

4's-5’s 

2.02E+01 

3.14E-06 
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Table  2. 


Ray 


eigh  Damping  Results  for  the  Vertical  Direction  by  Area  Group 


Area 

No. 

Deck 

Frame 

Athwartship 
Position 
p  =  port 
s  =  starboard 

Alpha 

Beta 

Number  of 
modes 
used  in 
curvefitting 
a  and  P 

1 

1 

60 

2'p 

2.12E+01 

3.54E-06 

12 

3 

1 

56 

4'p-2's 

1.96E+01 

4.22E-07 

33 

4 

1 

100 

15's 

1.95E+01 

4.08E-06 

26 

5 

1 

110 

CL 

1.97E+01 

1.62E-06 

55 

6 

1 

126.5 

CL-5's 

2.16E+01 

2.31E-07 

12 

7 

1 

142 

20’p-20’s 

2.09E+01 

2.11E-07 

104 

8 

1 

135 

17fs 

2.05E+01 

7.79E-06 

26 

9 

1 

159 

6’p-4's 

2.04E+01 

1.07E-06 

28 

10 

1 

171.5 

CL-15's 

2.03E+01 

1.16E-06 

41 

30 


Area 

No. 

Deck 

Frame 

Athwartship 
Position 
p  =  port 
s  =  starboard 

Alpha 

Beta 

Number  of 
modes 
used  in 
curvefitting 
a  and  P 

11 

1 

37.5 

5’p-CL 

2.04E+01 

2.77E-07 

16 

12 

1 

258 

8’p 

2.41E+01 

3.18E-05 

5 

13 

1 

312.5 

6’p-19s 

1.67E+01 

5.24E-07 

79 

14 

1 

324 

7'p-6'p 

1.96E+01 

2.82E-06 

21 

15 

1 

328 

5’p-5’s 

2.00E+01 

4.19E-07 

44 

16 

1 

375 

25’p-20’p 

1.79E+01 

8.38E-07 

129 

17 

1 

218 

28’p-CL 

1.55E+01 

1.91E-06 

72 

18 

2 

86.5 

i  rp-i  i’s 

1.85E+01 

1.82E-06 

69 

19 

2.5 

155 

3 1  fp-3 1 f  s 

2.03E+01 

-7.33E-08 

42 

20 

2 

137 

22’p-CL 

1.92E+01 

7.08E-07 

77 

21 

2 

134 

15’s 

1.63E+01 

2.06E-06 

22 

22 

2 

161 

12,p-3,s 

1.67E+01 

3.18E-06 

115 

23 

2 

162 

9's 

2.06E+01 

1.52E-06 

39 

24 

2 

306 

4’p-7’s 

1.87E+01 

1.53E-07 

19 

25 

2 

35 

3’p-ll's 

2.23E+01 

8.53E-07 

34 

26 

2 

452 

6’p 

1.94E+01 

3.00E-07 

28 

27 

2 

458.5 

12fs-22fs 

1.84E+01 

2.20E-06 

46 

28 

2 

70 

6’s-lO’s 

1.75E+01 

1.48E-06 

31 

29 

2 

75 

12'p-3’p 

2.11E+01 

2.17E-06 

19 

30 

2.5 

355.5 

10’p-l  l’s 

1.88E+01 

8.61E-07 

38 

31 

2.5 

398 

CL-8's 

2.07E+01 

1.51E-06 

67 

32 

3 

192.5 

6’p-23’s 

1.70E+01 

3.62E-06 

55 

33 

3 

277 

20'p-14's 

2.00E+01 

9.57E-07 

99 

34 

3 

119.5 

4'p-6's 

1.90E+01 

1.02E-06 

19 

35 

3 

129 

6’p-CL 

1.96E+01 

2.51E-06 

40 

36 

3 

149 

19'p-16’p 

2.16E+01 

3.84E-06 

5 

37 

3 

168 

12's 

2.22E+01 

1.73E-06 

19 

37 

3 

151 

36'p-36's 

1.63E+01 

9.65E-07 

54 

38 

3 

318-320 

2’p-15’s 

2.22E+01 

2.72E-07 

65 

40 

3 

331.5 

3's-CL 

2.07E+01 

1.86E-06 

47 

41 

3 

397.5 

13'p-CL 

1.50E+01 

2.77E-06 

32 

42 

3.5 

150 

19'p-6’s 

1.81E+01 

1.63E-06 

49 

43 

3.5 

156 

8’p-2fs 

2.11E+01 

1.50E-06 

49 

44 

3.5 

147 

16’p-15’0 

1.80E+01 

2.82E-06 

59 

46 

4 

150 

15's 

2.23E+01 

1.49E-06 

31 

47 

4 

156 

3’p-CL 

2.12E+01 

1.27E-07 

14 

48 

4 

174 

1  l’p-1  l’s 

1.90E+01 

2.65E-07 

25 

49 

4 

232.5 

6'p-CL 

1.84E+01 

1.01E-06 

30 

50 

4 

292 

CL 

1.66E+01 

2.13E-06 

16 

51 

4 

78 

2'P 

1.77E+01 

3.45E-06 

7 

31 


Area 

No. 

Deck 

Frame 

Athwartship 
Position 
p  =  port 
s  =  starboard 

Alpha 

Beta 

Number  of 
modes 
used  in 
curvefitting 
a  and  P 

52 

4 

96 

CL 

2.00E+01 

2.52E-05 

12 

53 

5 

280 

CL 

2.04E+01 

-5.41E-08 

42 

54 

5 

322.5 

17’p-15's 

1.91E+01 

3.35E-06 

26 

55 

6 

150 

CL 

2.25E+01 

3.10E-07 

17 

56 

6 

160 

CL 

2.03E+01 

5.66E-05 

12 

57 

7 

150 

CL 

2.35E+01 

5.21E-06 

7 

60 

8 

173.5 

32’p-32's 

1.01E+01 

6.16E-06 

10 

61 

9 

168 

CL 

1.90E+01 

9.63E-06 

10 

62 

9 

171.5 

25’p-25's 

2.09E+01 

4.86E-07 

89 

63 

10 

174 

3's-CL 

2.29E+01 

-1.24E-07 

16 

65 

11 

181 

CL 

2.01E+01 

4.57E-06 

15 

66 

HOLD 

433 

2'P 

1.29E+01 

1.24E-05 

10 

67 

MD 

281.5 

4's-5's 

2.02E+01 

2.84E-05 

4 

Table  3  shows  the  weighted  mean  value  of  a  for  each  of  the  directions  of  motion. 
The  weighting  factor,  shown  in  the  rightmost  column  of  Table  1,  is  defined  as  the  number 
of  modes  used  in  the  curve-fitting  a  and  (1.  Thus,  it  can  be  concluded  that  weighted  mean 
values  of  a  are  similar  in  both  directions. 


Table  4  shows  the  weighted  mean  value  of  p  for  each  direction  of  motion,  the 
weighting  factor  is  the  rightmost  column  of  Table  2,  the  number  of  modes  used  in  curve¬ 
fitting  a  and  p.  It  can  be  concluded  that  weighted  mean  values  in  athwartship  direction 
are  slightly  larger  than  those  of  the  vertical  direction.  Equation  68  is  used  to  calculate  the 
mean  a  coefficient, 


a 


mean 


Xa<N! 

i= 1 _ 

M 


(68) 
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where  the  variable  M  stand  for  the  of  areas  to  be  considered  . 


Table  3.  Weighted  Mean  of  a 


Athwartship  Direction 

Vertical  Direction 

18.4 

19.2 

Table  4.  Weighted  Mean  of  [3 
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2.82E-06 
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Figure  12.  Rayleigh  Damping  Coefficient  a  for  Athwartship  Direction  on  Deck  1 
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Figure  13.  Rayleigh  Damping  Coefficient  a  for  Vertical  Direction  on  Deck  1 
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V.  EFFECTS  OF  DAMPING  TO  SHIP  SHOCK  RESPONSES 


In  the  1994,  the  USS  JOHN  PAUL  JONES  (DDG  53)  was  chosen  as  the 
representative  ship  of  the  DDG  51  Arleigh  Burke  Class  Destroyer  and  subsequently 
subjected  to  a  series  of  shock  trials.  Some  seven  years  later  in  the  summer  of  2001, 
similar  ship  shock  trials  were  conducted  on  the  USS  WINSTON  S.  CHURCHILL  (DDG 
81).  This  latter  series  of  live  fire  tests  was  performed  on  the  31st  ship  in  the  same  class 
due  to  the  significant  design  changes  incorporated  into  the  Flight  IIA  version  of  this  ship. 
Some  of  the  significant  changes  that  were  found  in  DDG  8 1  included  an  extension  in  the 
ship’s  overall  length  and  the  additional  of  two  helicopter  hangers. 


Starting  with  a  highly  complex  finite  element  model  of  the  ship  and  the 
surrounding  fluid  mesh,  shown  in  Figure  14,  the  shock  response  velocity  was  calculated 
for  various  locations  throughout  the  ship  using  the  modeling  and  simulation  process 
outlined  in  Figure  15  [4], 


Figure  14.  DDG  81  Coupled  Fluid-Structure  Model 
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The  results  of  this  process  were  compared  with  the  actual  ship  shock  trial  sensor 
data  obtained  during  the  2001  Live  Fire  Testing  and  Evaluation. 
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Figure  15.  Model  Generation  and  Simulation  Flow  Chart 


Additionally  the  shock  velocity  response  plots  obtained  from  the  aforementioned 
process  were  compared  against  the  earlier  conducted  simulations  that  used  the  Rayleigh 
Damping  Coefficients  previously  used  in  the  DDG  53  modeling  and  simulation  effort, 
which  occurred  during  the  mid  1990’s  at  the  Naval  Postgraduate  School.  The  values  for 
these  coefficients  are  listed  in  Table  5. 
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Table  5.  Rayleigh  Damping  Coefficient  for  4%  &  8%  Proportional  Damping 


Damping  Value 

a 

P 

4% 

2.64 

4.99E-5 

8% 

4.93 

9.89E-5 

In  this  case  the  Rayleigh  parameters  (a,  P)  were  based  on  evaluation  of  the 
damping  values  at  two  given  points,  5  Hz  and  250  Hz,  which  cover  the  range  of  data 
which  was  required  for  subsequent  comparisons. 


RAYLEIGH  DAMPING 


Figure  16.  Rayleigh  Damping  Values  (in  Linear  Scale) 
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Figure  17.  Rayleigh  Damping  Values  (in  Logarithmic  Scale) 


The  following  series  of  velocity  response  plots  compares  the  Rayleigh  damping 
coefficients,  a  and  (3,  presented  in  Tables  3  and  4  with  coefficients  that  were  used  in 
previous  studies  conducted  on  the  DDG  53  and  DDG  81,  which  appear  in  Table  5. 


Russell’s  Error  Factor  [5-7]  was  chosen  as  a  means  of  comparing  the  velocity 
response  data  against  the  actual  ship  shock  trial  data.  For  the  purpose  of  this  study,  an 
established  set  of  acceptance  criteria  was  taken  from  the  work  accomplished  in  2003  on 
the  DDG  81  Ship  Shock  Simulation  [4],  These  values  are  presented  in  Table  6. 


Table  6.  Russell’s  Correlation  Acceptance  Criteria 


RC  <0.15 

Excellent 

0.15  <RC<  0.28 

Good 

RC  >  0.28 

Poor 
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Vertical  Velocity  (ft/s)  Vertical  Velocity  (ft/s) 


DDG-81  SHOT  2 

Grid  1 42489-vz  (V2002V) 


Time  (msec) 

Ship  Shock  Trial  Measured  Data 

NTS  Simulation  (alpha  =  19.2,  beta  =  2. 09  R- Ob) 
NTS  Simulation  (4%  Damping) 


Figure  18.  Sample  Vertical  Velocity  Response:  Deck  Sensor 


DDG-81  SHOT  2 


Grid  210894  -vz  (V2008VI) 


Time  (msec) 

Ship  Shock  Trial  Measured  Data 

NTS  Simulation  (alpha  =  19.2,  beta  =  2.09R-06) 
NTS  Simulation  (4%  Damping) 


Figure  19.  Sample  Vertical  Velocity  Response:  Deck  Sensor 
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Vertical  Velocity  (ft/s)  Vertical  Velocity  (ft/s) 


DDG-8 1  SHOT  2 

Grid  330769-vz  (V2035V) 


Time  (msec) 

Ship  Shock  Trial  Measured  Data 

NPS  Simulation  (alpha  =  19.2,  beta  =  2. 09  R- 06) 
NTS  Simulation  (4%  Damping) 


Figure  20.  Sample  Vertical  Velocity  Response:  Keel  Sensor 


DDG-8 1  SHOT  2 

Grid  222436-vz  (V2125V) 


Time  (msec) 

Ship  Shock  Trial 

NTS  Simulation  (alpha  =  19.2,  beta  =  2.09R-06) 
NPS  Simulation  (4%  Damping) 


Figure  21.  Sample  Vertical  Velocity  Response:  Bulkhead  Sensor 
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As  the  velocity  response  plot  comparisons  in  Figure  18  through  Figure  21  show, 
there  is  a  much  better  correlation  between  the  NPS  damping  values  and  the  ship  shock 
trial  data,  than  with  the  fixed  4%  damping.  For  the  sensors  examined  of  which  the 
approximate  location  of  each  is  indicated  on  the  time  history  plots  by  a  red  dot,  a 
Russell’s  Comprehensive  (RC)  error  correlation  factor  was  computed.  The  mean  RC  for 
the  4%  Damping  cases  was  0.25  while  in  comparison  when  the  new  NPS  damping  values 
from  Table  3  and  Table  4  were  used,  the  mean  RC  value  was  only  0.18.  Recalling  that 
by  Russell’s  correlation  criteria,  a  value  below  0.15  is  considered  an  excellent 
correlation,  the  simulations  using  the  new  NPS  damping  values  consistently  show  better 
overall  correlation  and  an  average  reduction  of  approximately  25%  in  deviation  from  the 
recorded  ship  shock  trial  data  versus  those  using  the  fixed  4%  damping.  Table  7 
illustrates  a  sampling  of  the  supporting  data. 


Table  7.  Comparison  of  Russell’s  Error  Factor  for  DDG  8 1  Shot  2  (vertical  direction) 


Ship  Shock  Simulation  with  Shot  2  Geometry,  Dense  Mesh  and  738  in  Cavitation  Depth 


Shock  Trial  Data  vs.  4%  Damping 

Shock  Trial  Data  vs.  NPS  Damping 

Sensor 

Node 

Mounting 

Location  (in)* 

SHIP  SHOCK  TRIAL  DATA  (<250HZ) 

SHIP  SHOCK  TRIAL  DATA  (<250HZ) 

Type 

X 

y 

z 

RM 

RP 

RC 

RM 

RP 

RC 

■EMBlil 

Deck 

24 

85 

0.1974 

0.2715 

0.2975 

0.0679 

0.2175 

0.2019 

Deck 

176 

171 

0.1207 

0.2689 

0.258 

0.1200 

0.1932 

0.2016 

■cgTOEEl 

Keel 

1152 

135 

193 

0.1643 

0.1849 

0.2192 

0.0009 

0.1692 

0.15 

Bulkhead 

375 

390 

0.1651 

0.1936 

0.2255 

0.0214 

0.1914 

0.1707 

Russell  Error  Correlation 

Sum(E(X)) 

0.64750 

0.91890 

1.00020 

0.21020 

0.77130 

0.72420 

>0.28 

Poor 

Sum(E(XA2)) 

0.10779 

0.21769 

0.25397 

0.01947 

0.14990 

0.13304 

<0.15 

Excellent 

Mean 

0.16188 

0.22973 

0.25005 

0.05255 

0.19283 

0.18105 

Standard  Deviation 

0.03148 

0.04688 

0.03591 

0.05299 

0.01974 

0.02535 

*  Referenced  to  the  G&C  NASTRAN  Model  coordinate  origin  located  at  the  stern.  In  the  Y-direction,  port  is  postive  from  the  centerline. 


Figure  22  is  a  graphical  representation  of  the  data  presented  in  Table  7.  Notice 
that  the  Russell’s  Comprehensive  correlation  factor  for  the  simulations  using  the  NPS 
damping  values  are  all  in  the  excellent  or  highly  acceptable  range,  while  the  results  from 
the  simulations  performed  using  the  fixed  4%  damping  values  are  only  marginally 
acceptable  or  fall  outside  of  the  acceptable  range  all  together.  Note  as  well  that  there  is 
considerable  improvement  in  the  accuracy  of  the  magnitude  component  of  the  Russell’s 
correlation  in  the  simulations  using  the  NPS  damping  values,  as  demonstrated  by  the 
grouping  of  points  nearer  the  ordinate. 
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DDG  81  SHOT  2 


- 0.15  - 0.28  ♦  NPS  Damping  Values  •  4%  Damping 


Figure  22.  Russell’s  Error  Factor  for  Selected  Sensors  of  DDG  81  Shot  2 

In  this  small  but  representative  sampling  of  data  points  from  various  locations 
throughout  the  ship,  comparison  of  the  Russell’s  error  correlation  shows  that  using  the 
NPS  damping  values  tends  to  improve  the  accuracy  of  the  simulation  from  20%  to  30%. 
These  results  are  presented  in  Table  8. 


Table  8.  Relative  %  Change  in  RC  for  NPS  Damping  versus  4%  Damping  Case 


Sensor 

Node 

Ship 

Compartment 

Location 

Percent 

Relavtive 

Change 

V2008VI 

210894 

4th  Deck 

32.13% 

V2002V 

142489 

4th  Deck 

21.86% 

V2035V 

330769 

3rd  Deck 

31.57% 

V2125V 

222436 

1st  Deck 

24.30% 

Average  Improvement  in  Correlation 

27.47% 
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VI.  CONCLUSIONS 


Rayleigh  damping  representation  for  modeling  ship  system  damping  has  been 
investigated  based  on  the  ship  shock  trial  data.  Based  on  the  results  of  studies,  a  set  of 
Rayleigh  damping  parameters  (a  and  (3)  are  recommended  for  ship  shock  response 
predictions.  The  results  of  investigation  also  indicate  that  (i)  the  system  damping  is 
largely  affected  by  mass  driven  (a[M]),  and  (ii)  the  damping  decreases  as  frequency 
increases  as  we  commonly  understood. 
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APPENDIX  A.  FIGURES  OF  THE  RESULTS  OF  RAYLEIGH  DAMPING 
CURVE-FITTING 

A.  Results  in  the  Athwartship  Direction 


Figure  23.  Modal  Damping  Ratio  at  Area  1,  Athwartship  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  24.  Modal  Damping  Ratio  at  Area  1,  Athwartship  Direction  (Modified) 


Frequency  (Hz) 


Figure  25.  Modal  Damping  Ratio  at  Area  3,  Athwartship  Direction  (Original) 
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Figure  26.  Modal  Damping  Ratio  at  Area  3,  Athwartship  Direction  (Modified) 


Figure  27.  Modal  Damping  Ratio  at  Area  6,  Athwartship  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  28.  Modal  Damping  Ratio  at  Area  6,  Athwartship  Direction  (Modified) 
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Figure  29.  Modal  Damping  Ratio  at  Area  7,  Athwartship  Direction  (Original) 
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Figure  30. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  7,  Athwartship  Direction  (Modified) 


Figure  3 1 .  Modal  Damping  Ratio  at  Area  8,  Athwartship  Direction  (Original) 
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Damping  Ratio  Damping  Ratio 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  32.  Modal  Damping  Ratio  at  Area  8,  Athwartship  Direction  (Modified) 


Frequency  (Hz) 


Figure  33.  Modal  Damping  Ratio  at  Area  10,  Athwartship  Direction  (Original) 
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Figure  34.  Modal  Damping  Ratio  at  Area  10,  Athwartship  Direction  (Modified) 
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Figure  35.  Modal  Damping  Ratio  at  Area  11,  Athwartship  Direction  (Original) 
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Figure  36. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  11,  Athwartship  Direction  (Modified) 


Figure  37.  Modal  Damping  Ratio  at  Area  12,  Athwartship  Direction  (Original) 
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Figure  38.  Modal  Damping  Ratio  at  Area  12,  Athwartship  Direction  (Modified) 


Figure  39.  Modal  Damping  Ratio  at  Area  15,  Athwartship  Direction  (Original) 


E 

a  0-&1 


Measured 

CwrveFil 


Frequency  (Hz) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  40.  Modal  Damping  Ratio  at  Area  15,  Athwartship  Direction  (Modified) 
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Figure  41.  Modal  Damping  Ratio  at  Area  16,  Athwartship  Direction  (Original) 
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(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  42.  Modal  Damping  Ratio  at  Area  16,  Athwartship  Direction  (Modified) 


Figure  43.  Modal  Damping  Ratio  at  Area  18,  Athwartship  Direction  (Original) 
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Figure  44.  Modal  Damping  Ratio  at  Area  18,  Athwartship  Direction  (Modified) 


Figure  45.  Modal  Damping  Ratio  at  Area  19,  Athwartship  Direction  (Original) 


Figure  46. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  19,  Athwartship  Direction  (Modified) 
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Figure  47.  Modal  Damping  Ratio  at  Area  20,  Athwartship  Direction  (Original) 
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Figure  48.  Modal  Damping  Ratio  at  Area  20,  Athwartship  Direction  (Modified) 


Figure  49.  Modal  Damping  Ratio  at  Area  28,  Athwartship  Direction  (Original) 
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(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  50.  Modal  Damping  Ratio  at  Area  28,  Athwartship  Direction  (Modified) 
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Figure  5 1 .  Modal  Damping  Ratio  at  Area  30,  Athwartship  Direction  (Original) 


Figure  52. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  30,  Athwartship  Direction  (Modified) 
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Figure  53.  Modal  Damping  Ratio  at  Area  3 1 ,  Athwartship  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  54.  Modal  Damping  Ratio  at  Area  31,  Athwartship  Direction  (Modified) 


Figure  55.  Modal  Damping  Ratio  at  Area  32,  Athwartship  Direction  (Original) 
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Figure  56.  Modal  Damping  Ratio  at  Area  32,  Athwartship  Direction  (Modified) 


Figure  57.  Modal  Damping  Ratio  at  Area  33,  Athwartship  Direction  (Original) 


50  IDO  150  200  250 

Frequency  (Hz) 


E 

O  0-&1 


11  " 

■ 

■ 

*  ^  . 

* 

■ 

■  »  V  *  ■ 

"  -  .V  * 1  - 

p"  *  i  > 

"  *!*■  STV.' 

m 

mm 

** 

■ 

■ 

■  Measured 

CucveFii 

Frequency  (Hz) 


Figure  58. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  33,  Athwartship  Direction  (Modified) 


57 


1.0 


<*  0-5 

t  •* 


Measured 

Gurvefil 


100  ISO 

Frequency  {Hi) 


Figure  59.  Modal  Damping  Ratio  at  Area  35,  Athwartship  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  60.  Modal  Damping  Ratio  at  Area  35,  Athwartship  Direction  (Modified) 
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Figure  61.  Modal  Damping  Ratio  at  Area  37,  Athwartship  Direction  (Original) 
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Frequency  (Hz) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  62.  Modal  Damping  Ratio  at  Area  37,  Athwartship  Direction  (Modified) 


Figure  63.  Modal  Damping  Ratio  at  Area  40,  Athwartship  Direction  (Original) 


Figure  64. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  40,  Athwartship  Direction  (Modified) 
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Figure  65.  Modal  Damping  Ratio  at  Area  41,  Athwartship  Direction  (Original) 
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Figure  66. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  41,  Athwartship  Direction  (Modified) 


Figure  67.  Modal  Damping  Ratio  at  Area  42,  Athwartship  Direction  (Original) 
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(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  68.  Modal  Damping  Ratio  at  Area  42,  Athwartship  Direction  (Modified) 
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Figure  69.  Modal  Damping  Ratio  at  Area  43,  Athwartship  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  70.  Modal  Damping  Ratio  at  Area  43,  Athwartship  Direction  (Modified) 
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Figure  7 1 .  Modal  Damping  Ratio  at  Area  47,  Athwartship  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  72.  Modal  Damping  Ratio  at  Area  47,  Athwartship  Direction  (Modified) 


Figure  73.  Modal  Damping  Ratio  at  Area  48,  Athwartship  Direction  (Original) 
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(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  74.  Modal  Damping  Ratio  at  Area  48,  Athwartship  Direction  (Modified) 


Figure  75.  Modal  Damping  Ratio  at  Area  49,  Athwartship  Direction  (Original) 
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Figure  76.  Modal  Damping  Ratio  at  Area  49,  Athwartship  Direction  (Modified) 
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Figure  77.  Modal  Damping  Ratio  at  Area  50,  Athwartship  Direction  (Original) 
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Figure  78.  Modal  Damping  Ratio  at  Area  50,  Athwartship  Direction  (Modified) 


Figure  79.  Modal  Damping  Ratio  at  Area  53,  Athwartship  Direction  (Original) 
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(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  80.  Modal  Damping  Ratio  at  Area  53,  Athwartship  Direction  (Modified) 


Figure  81.  Modal  Damping  Ratio  at  Area  54,  Athwartship  Direction  (Original) 
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Figure  82. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  54,  Athwartship  Direction  (Modified) 
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Figure  83.  Modal  Damping  Ratio  at  Area  55,  Athwartship  Direction  (Original) 
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Figure  84.  Modal  Damping  Ratio  at  Area  55,  Athwartship  Direction  (Modified) 


Figure  85.  Modal  Damping  Ratio  at  Area  56,  Athwartship  Direction  (Original) 
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Figure  86.  Modal  Damping  Ratio  at  Area  56,  Athwartship  Direction  (Modified) 


Figure  87.  Modal  Damping  Ratio  at  Area  57,  Athwartship  Direction  (Original) 
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Figure  88.  Modal  Damping  Ratio  at  Area  57,  Athwartship  Direction  (Modified) 
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Figure  89.  Modal  Damping  Ratio  at  Area  58,  Athwartship  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  90.  Modal  Damping  Ratio  at  Area  58,  Athwartship  Direction  (Modified) 


Figure  91.  Modal  Damping  Ratio  at  Area  59,  Athwartship  Direction  (Original) 
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Figure  92.  Modal  Damping  Ratio  at  Area  59,  Athwartship  Direction  (Modified) 


Figure  93.  Modal  Damping  Ratio  at  Area  60,  Athwartship  Direction  (Original) 


Figure  94. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  60,  Athwartship  Direction  (Modified) 
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Figure  95.  Modal  Damping  Ratio  at  Area  61,  Athwartship  Direction  (Original) 
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Figure  96.  Modal  Damping  Ratio  at  Area  61,  Athwartship  Direction  (Modified) 
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Figure  97.  Modal  Damping  Ratio  at  Area  62,  Athwartship  Direction  (Original) 
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(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  98.  Modal  Damping  Ratio  at  Area  62,  Athwartship  Direction  (Modified) 


Figure  99.  Modal  Damping  Ratio  at  Area  64,  Athwartship  Direction  (Original) 
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Figure  100.  Modal  Damping  Ratio  at  Area  64,  Athwartship  Direction  (Modified) 
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Figure  101.  Modal  Damping  Ratio  at  Area  65,  Athwartship  Direction  (Original) 
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Figure  102.  Modal  Damping  Ratio  at  Area  65,  Athwartship  Direction  (Modified) 


Figure  103.  Modal  Damping  Ratio  at  Area  67,  Athwartship  Direction  (Original) 
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Figure  104. 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Modal  Damping  Ratio  at  Area  67,  Athwartship  Direction  (Modified) 
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B.  Results  in  the  Vertical  Direction 


Figure  105.  Modal  Damping  Ratio  at  Area  1,  Vertical  Direction  (Original) 


(in  Linear  Scale) 


(in  Logarithmic  Scale) 


Figure  106.  Modal  Damping  Ratio  at  Area  1,  Vertical  Direction  (Modified) 


Figure  107.  Modal  Damping  Ratio  at  Area  3,  Vertical  Direction  (Original) 
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Figure  108.  Modal  Damping  Ratio  at  Area  3,  Vertical  Direction  (Modified) 


Figure  109.  Modal  Damping  Ratio  at  Area  4,  Vertical  Direction  (Original) 
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Figure  1 10.  Modal  Damping  Ratio  at  Area  4,  Vertical  Direction  (Modified) 
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Figure  111.  Modal  Damping  Ratio  at  Area  5,  Vertical  Direction  (Original) 
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Figure  112.  Modal  Damping  Ratio  at  Area  5,  Vertical  Direction  (Modified) 
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Figure  1 13.  Modal  Damping  Ratio  at  Area  6,  Vertical  Direction  (Original) 
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Figure  1 14.  Modal  Damping  Ratio  at  Area  6,  Vertical  Direction  (Modified) 
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Figure  115.  Modal  Damping  Ratio  at  Area  7,  Vertical  Direction  (Original) 
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Figure  1 16.  Modal  Damping  Ratio  at  Area  7,  Vertical  Direction  (Modified) 
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Figure  117.  Modal  Damping  Ratio  at  Area  8,  Vertical  Direction  (Original) 
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Figure  118.  Modal  Damping  Ratio  at  Area  8,  Vertical  Direction  (Modified) 
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Figure  1 19.  Modal  Damping  Ratio  at  Area  9,  Vertical  Direction  (Original) 
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Figure  120.  Modal  Damping  Ratio  at  Area  9,  Vertical  Direction  (Modified) 


Figure  121.  Modal  Damping  Ratio  at  Area  10,  Vertical  Direction  (Original) 


J 

1 

I 

Tht 

mm 

-  •  - 

■  ,  ■ 

h  m 

D  "  ■  « 

"  ^  "j 
pl|  « 

V 

v 

■ 

■  Measured 

CurveFit 

10  too 


Frequency  (Hz) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  122.  Modal  Damping  Ratio  at  Area  10,  Vertical  Direction  (Modified) 
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Figure  123.  Modal  Damping  Ratio  at  Area  11,  Vertical  Direction  (Original) 
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Figure  124.  Modal  Damping  Ratio  at  Area  11,  Vertical  Direction  (Modified) 


Figure  125.  Modal  Damping  Ratio  at  Area  12,  Vertical  Direction  (Original) 
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Figure  126.  Modal  Damping  Ratio  at  Area  12,  Vertical  Direction  (Modified) 
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Figure  127.  Modal  Damping  Ratio  at  Area  13,  Vertical  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  128.  Modal  Damping  Ratio  at  Area  13,  Vertical  Direction  (Modified) 
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Figure  129.  Modal  Damping  Ratio  at  Area  14,  Vertical  Direction  (Original) 
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Figure  130.  Modal  Damping  Ratio  at  Area  14,  Vertical  Direction  (Modified) 


Figure  131.  Modal  Damping  Ratio  at  Area  15,  Vertical  Direction  (Original) 
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Figure  132.  Modal  Damping  Ratio  at  Area  15,  Vertical  Direction  (Modified) 
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Figure  133.  Modal  Damping  Ratio  at  Area  16,  Vertical  Direction  (Original) 
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Figure  134.  Modal  Damping  Ratio  at  Area  16,  Vertical  Direction  (Modified) 
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Figure  135.  Modal  Damping  Ratio  at  Area  17,  Vertical  Direction  (Original) 
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Figure  136.  Modal  Damping  Ratio  at  Area  17,  Vertical  Direction  (Modified) 
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Figure  137.  Modal  Damping  Ratio  at  Area  18,  Vertical  Direction  (Original) 
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Figure  138.  Modal  Damping  Ratio  at  Area  18,  Vertical  Direction  (Modified) 
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Figure  139.  Modal  Damping  Ratio  at  Area  19,  Vertical  Direction  (Original) 
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Figure  140.  Modal  Damping  Ratio  at  Area  19,  Vertical  Direction  (Modified) 
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Figure  141.  Modal  Damping  Ratio  at  Area  20,  Vertical  Direction  (Original) 
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Figure  142.  Modal  Damping  Ratio  at  Area  20,  Vertical  Direction  (Modified) 


Figure  143.  Modal  Damping  Ratio  at  Area  21,  Vertical  Direction  (Original) 
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Figure  144.  Modal  Damping  Ratio  at  Area  21,  Vertical  Direction  (Modified) 
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Figure  145.  Modal  Damping  Ratio  at  Area  22,  Vertical  Direction  (Original) 
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Figure  146.  Modal  Damping  Ratio  at  Area  22,  Vertical  Direction  (Modified) 
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Figure  147.  Modal  Damping  Ratio  at  Area  23,  Vertical  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  148.  Modal  Damping  Ratio  at  Area  23,  Vertical  Direction  (Modified) 


Figure  149.  Modal  Damping  Ratio  at  Area  24,  Vertical  Direction  (Original) 
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Figure  150.  Modal  Damping  Ratio  at  Area  24,  Vertical  Direction  (Modified) 


Figure  151.  Modal  Damping  Ratio  at  Area  25,  Vertical  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  152.  Modal  Damping  Ratio  at  Area  25,  Vertical  Direction  (Modified) 
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Figure  153.  Modal  Damping  Ratio  at  Area  26,  Vertical  Direction  (Original) 
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Figure  154.  Modal  Damping  Ratio  at  Area  26,  Vertical  Direction  (Modified) 
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Figure  155.  Modal  Damping  Ratio  at  Area  27,  Vertical  Direction  (Original) 
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Figure  156.  Modal  Damping  Ratio  at  Area  27,  Vertical  Direction  (Modified) 


Figure  157.  Modal  Damping  Ratio  at  Area  28,  Vertical  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  158.  Modal  Damping  Ratio  at  Area  28,  Vertical  Direction  (Modified) 
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Figure  159.  Modal  Damping  Ratio  at  Area  29,  Vertical  Direction  (Original) 
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Figure  160.  Modal  Damping  Ratio  at  Area  29,  Vertical  Direction  (Modified) 
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Figure  161 .  Modal  Damping  Ratio  at  Area  30,  Vertical  Direction  (Original) 


92 


■  Measured 

CurveFit 

■ 

.1 

• . 

_  ■ 

■  ■ 

-  _ . 

50  100  150  200  250 


Frequency  (Hz) 


Frequency  (Hz) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  162.  Modal  Damping  Ratio  at  Area  30,  Vertical  Direction  (Modified) 


Figure  163.  Modal  Damping  Ratio  at  Area  31,  Vertical  Direction  (Original) 
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Figure  164.  Modal  Damping  Ratio  at  Area  31,  Vertical  Direction  (Modified) 
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Figure  165.  Modal  Damping  Ratio  at  Area  32,  Vertical  Direction  (Original) 
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Figure  166.  Modal  Damping  Ratio  at  Area  32,  Vertical  Direction  (Modified) 
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Figure  167.  Modal  Damping  Ratio  at  Area  33,  Vertical  Direction  (Original) 


94 


Damping  Ratio  Damping  Ratio 


Frequency  (Hz) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  168.  Modal  Damping  Ratio  at  Area  33,  Vertical  Direction  (Modified) 


Figure  169.  Modal  Damping  Ratio  at  Area  34,  Vertical  Direction  (Original) 
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Figure  170.  Modal  Damping  Ratio  at  Area  34,  Vertical  Direction  (Modified) 
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Figure  171.  Modal  Damping  Ratio  at  Area  35,  Vertical  Direction  (Original) 
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Figure  172.  Modal  Damping  Ratio  at  Area  35,  Vertical  Direction  (Modified) 
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Figure  173.  Modal  Damping  Ratio  at  Area  36,  Vertical  Direction  (Original) 
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Figure  174.  Modal  Damping  Ratio  at  Area  36,  Vertical  Direction  (Modified) 


Figure  175.  Modal  Damping  Ratio  at  Area  37,  Vertical  Direction  (Original) 
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Figure  176.  Modal  Damping  Ratio  at  Area  37,  Vertical  Direction  (Modified) 
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Figure  177.  Modal  Damping  Ratio  at  Area  38,  Vertical  Direction  (Original) 
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Figure  178.  Modal  Damping  Ratio  at  Area  38,  Vertical  Direction  (Modified) 


Figure  179.  Modal  Damping  Ratio  at  Area  39,  Vertical  Direction  (Original) 
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Figure  180.  Modal  Damping  Ratio  at  Area  39,  Vertical  Direction  (Modified) 


Figure  181.  Modal  Damping  Ratio  at  Area  40,  Vertical  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  182.  Modal  Damping  Ratio  at  Area  40,  Vertical  Direction  (Modified) 
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Figure  183.  Modal  Damping  Ratio  at  Area  41,  Vertical  Direction  (Original) 
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Figure  184.  Modal  Damping  Ratio  at  Area  41,  Vertical  Direction  (Modified) 
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Figure  185.  Modal  Damping  Ratio  at  Area  42,  Vertical  Direction  (Original) 
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Figure  186.  Modal  Damping  Ratio  at  Area  42,  Vertical  Direction  (Modified) 


Figure  187.  Modal  Damping  Ratio  at  Area  43,  Vertical  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  188.  Modal  Damping  Ratio  at  Area  43,  Vertical  Direction  (Modified) 
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Figure  189.  Modal  Damping  Ratio  at  Area  44,  Vertical  Direction  (Original) 
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(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  190.  Modal  Damping  Ratio  at  Area  44,  Vertical  Direction  (Modified) 


Figure  191 .  Modal  Damping  Ratio  at  Area  45,  Vertical  Direction  (Original) 
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Figure  192.  Modal  Damping  Ratio  at  Area  45,  Vertical  Direction  (Modified) 


Figure  193.  Modal  Damping  Ratio  at  Area  46,  Vertical  Direction  (Original) 
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Figure  194.  Modal  Damping  Ratio  at  Area  46,  Vertical  Direction  (Modified) 
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Figure  195.  Modal  Damping  Ratio  at  Area  47,  Vertical  Direction  (Original) 
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Figure  196.  Modal  Damping  Ratio  at  Area  47,  Vertical  Direction  (Modified) 


Figure  197.  Modal  Damping  Ratio  at  Area  48,  Vertical  Direction  (Original) 
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Figure  198.  Modal  Damping  Ratio  at  Area  48,  Vertical  Direction  (Modified) 


Figure  199.  Modal  Damping  Ratio  at  Area  49,  Vertical  Direction  (Original) 
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Figure  200.  Modal  Damping  Ratio  at  Area  49,  Vertical  Direction  (Modified) 
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Figure  201.  Modal  Damping  Ratio  at  Area  50,  Vertical  Direction  (Original) 
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Figure  202.  Modal  Damping  Ratio  at  Area  50,  Vertical  Direction  (Modified) 


Figure  203.  Modal  Damping  Ratio  at  Area  5 1 ,  Vertical  Direction  (Original) 
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Figure  204.  Modal  Damping  Ratio  at  Area  51,  Vertical  Direction  (Modified) 


Figure  205.  Modal  Damping  Ratio  at  Area  52,  Vertical  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  206.  Modal  Damping  Ratio  at  Area  52,  Vertical  Direction  (Modified) 
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Figure  207.  Modal  Damping  Ratio  at  Area  53,  Vertical  Direction  (Original) 
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(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  208.  Modal  Damping  Ratio  at  Area  53,  Vertical  Direction  (Modified) 
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Figure  209.  Modal  Damping  Ratio  at  Area  54,  Vertical  Direction  (Original) 
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Figure  210.  Modal  Damping  Ratio  at  Area  54,  Vertical  Direction  (Modified) 


Figure  211.  Modal  Damping  Ratio  at  Area  55,  Vertical  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  212.  Modal  Damping  Ratio  at  Area  55,  Vertical  Direction  (Modified) 
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Figure  213.  Modal  Damping  Ratio  at  Area  56,  Vertical  Direction  (Original) 
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(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  214.  Modal  Damping  Ratio  at  Area  56,  Vertical  Direction  (Modified) 


Figure  215.  Modal  Damping  Ratio  at  Area  57,  Vertical  Direction  (Original) 
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Figure  216.  Modal  Damping  Ratio  at  Area  57,  Vertical  Direction  (Modified) 


Figure  217.  Modal  Damping  Ratio  at  Area  59,  Vertical  Direction  (Original) 


(in  Linear  Scale)  (in  Logarithmic  Scale) 

Figure  218.  Modal  Damping  Ratio  at  Area  59,  Vertical  Direction  (Modified) 


111 


1.0 


09 
08 
07 

2  06 
5 

05 
04 
03 
0  2 
01 
00 

0  50  100  150  200  250 

Frequency  (Hz) 

Figure  219.  Modal  Damping  Ratio  at  Area  60,  Vertical  Direction  (Original) 
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Figure  220.  Modal  Damping  Ratio  at  Area  60,  Vertical  Direction  (Modified) 


Figure  22 1 .  Modal  Damping  Ratio  at  Area  6 1 ,  Vertical  Direction  (Original) 
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Figure  222.  Modal  Damping  Ratio  at  Area  61,  Vertical  Direction  (Modified) 


Figure  223.  Modal  Damping  Ratio  at  Area  62,  Vertical  Direction  (Original) 
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Figure  224.  Modal  Damping  Ratio  at  Area  62,  Vertical  Direction  (Modified) 
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Figure  225.  Modal  Damping  Ratio 
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Figure  226.  Modal  Damping  Ratio 
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Figure  227.  Modal  Damping  Ratio  at  Area  64,  Vertical  Direction  (Original) 
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Figure  228.  Modal  Damping  Ratio  at  Area  64,  Vertical  Direction  (Modified) 


Figure  229.  Modal  Damping  Ratio  at  Area  65,  Vertical  Direction  (Original) 
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Figure  230.  Modal  Damping  Ratio  at  Area  65,  Vertical  Direction  (Modified) 
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Figure  23 1 .  Modal  Damping  Ratio  at  Area  66,  Vertical  Direction  (Original) 
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Figure  232.  Modal  Damping  Ratio  at  Area  66,  Vertical  Direction  (Modified) 


Figure  233.  Modal  Damping  Ratio  at  Area  67,  Vertical  Direction  (Original) 
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Figure  234.  Modal  Damping  Ratio  at  Area  67,  Vertical  Direction  (Modified) 
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APPENDIX  B.  MODAL  PARAMETER  EXTRACTION  PROGRAM  LIST  IN 
TIME  DOMAIN 

This  program  uses  IMSL  libraries  that  are  included  in  Microsoft  Power  Station. 

use  msimslmd 

implicit  double  precision  (a-h,  o-z) 

dimension  tt  (4000) ,  w(4000),  hh  (4000, 200) ,  hh2n(4000),  betaq(201),  coeff(201),  & 
coefre (0 : 201) ,  coefim(0 : 201 ) ,  rootre(200),  rootim(200) 
dimension  hht (200, 4000) ,hhthh (200, 200) ,hht2n (200) 
double  complex  root  (200),  freg(6000),  wregm(6000) 
dimension  xr  (200) ,  wr(200),  wregr(4000),  wregi(4000) 
dimension  arn(200),  nonq(1000) 

double  complex  aaa (4000, 200) ,aaat (200, 4000) , ah (4000) ,aaatah (200) ,  & 
aaataaa (200, 200) ,  cmodal (201) , sr (201) , cw 
character*25  ifiles (1000) ,pfiles (1000) ,  gfiles (1000) , sensor (1000) ,  sensorid 
LOGICAL  COMPL 

pi=4 . 0d0*datan (1 . OdO) 
tfinal=2000 . OdO 

nf=773  !  no.  of  data  sets  774 

!  nq=170  !  total  mode  no.  to  be  considered 

open (uni t=l,  f ile=' modi fieddat. 1st 1 )  !  input  time  data  file  list 

do  1881  if  1=1.,  nf 

read(l,*)  sensor (ifl),  nonq(ifl) 

ifiles (ifl) =trim( sensor (ifl) ) // ’_2ms.dat ' 
pfiles (ifl) =trim  (sensor (ifl) ) //' _modal.dat  ’ 
gfiles (ifl) =trim (sensor (ifl) ) // '_regen.dat ’ 

1881  continue 

do  8888  if 1=4 67, nf 

mm=1500  !  data  points  to  be  used  for  calculation 
nq=nonq(ifl) 

write (*, *)  ifl, 1  — >  ',  ifiles (ifl) 
open(unit=3,  file=ifiles (ifl) ,  err=8888) 
open (unit=4,  file=pfiles (ifl) ) 
open (unit=6,  file=gfiles (ifl) ) 

i 

!  Velocity  Data  input,  time  is  in  miliseconds. 

i 

do  j=l,mm 

read (3, *)  tt(j),  w(j) 
if (tt (j) .ge.tfinal)  goto  1788 
enddo 

1788  mm=j-l 

i 

dt=  (tt  (mm)  -tt  (1) )  /dfloat  (mm-1)  /1000 .  OdO  !  time  interval  in  seconds 

!  write ( * , * )  dt 

i 

!  Calculate  [H] {beta}={h} 

i 

m2=mm-nq-l 
do  i=l,m2 
ii=i-l 

hh2n(i)=  -w(ii+nq+l) 
do  j=l,nq 

hh  (i,  j )  =w  (ii+j ) 
enddo 

enddo 
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hht  (l:nq,  1  :m2)  ^transpose  (hh  (1  :m2,  1  :nq) ) 

hhthh  (1  :nq,  1  :nq)  =matmul  (hht  (1  :nq,  1  :m2)  ,hh  (1  :m2, 1  :nq) ) 

hht2n(l:nq)  =matmul  (hht  (1  :nq,  1  :m2)  ,hh2n(l:m2) ) 

CALL  DLSARG  (nq,  hhthh,  200,  hht2n,  1,  betaq) 

i 

!  Solve  Polynomial  Equation,  and  calculate  {V} 

i 

betaq  (nq+1 )  =1 .  OdO 
compl= . FALSE . 

coef re ( 0 : nq) =betaq ( 1 : nq+1 ) 

CALL  BAUPOL  (coef  re,  coef  dm,  nq,  COMPL,  ROOTRE,  ROOTIM,  NUMIT) 
root  ( 1 :  nq)  =dcmplx  ( rootre  ( 1 :  nq) ,  root im  ( 1 :  nq) ) 

i 

!  Calculate  modal  frquencies  and  damping  ratios  from  {V} 

i 

do  i=l,nq 

sr (i) =cdlog (root (i) ) /dt 
wr  (i)  =cdabs  (sr  (i) ) 
xr(i)=  -dble  (sr  (i) ) /wr  (i) 

enddo 

i 

!  Calculate  modal  constants  {Ar} 

i 

cmodal=dcmplx  ( 0 .  OdO ,  0 .  OdO ) 
m2=mm-l 
do  i=l,m2 

ah  (i)  =w(i) 
do  j=l,nq 

aaa (i, j ) =root ( j ) **dfloat (1—1) 
enddo 

enddo 

aaat ( 1 : nq, 1 : m2 ) ^transpose (aaa ( 1 : m2 , 1 : nq) ) 

aaataaa  (1  :nq,  1  :nq)  =matmul  (aaat  (1  :nq,  1  :m2) ,  aaa  (1  :m2, 1  :nq) ) 

aaatah(linq)  =mtmul  (aaat  (1  :nq,  1  :m2)  ,ah(l:m2) ) 

CALL  DLSACG  (nq,  aaataaa,  200,  aaat ah,  1,  cmodal) 

i 

sensorid=sensor (ifl) 
do  i=l,nq 

cmodal (i) =cmodal (i) /sr (i) 
enddo 

i 

armax=0 .  OdO 
arn=0 . OdO 
do  i=l,nq 

if  (cdabs  (cmodal  (i) )  . gt .  armax)  armax=cdabs  (cmodal  (i) ) 

enddo 

arn  ( 1 :  nq)  =cdabs  ( cmodal  ( 1 :  nq) )  /armax 

i 

!  Calculate  Spectrum  using  modal  parameter  to  verify  results 

i 

fspan=l . OdO/dt/2 . OdO 
wspan=2 . OdO ^pi^f span 
dw=wspan/  (dfloat  (irm/2)  -1 .  OdO) 
f  reg^dcmplx  ( 0 .  OdO ,  0 .  OdO ) 
do  i=2,mm/2+l 

w=dfloat  (1—1 )  *dw 

cw=dcmplx  ( 0 .  OdO ,  w) 
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do  j=l,nq 

freg  (i)  =freg  (i)  +cmodal  ( j )  *cw/  (cw-sr  ( j ) ) 
enddo 

enddo 

do  i=2,  rrm/2 

ii=mm-i+2 

freg  (ii)  =dconjg  (freg  (i) ) 
enddo 

freg=freg*dw/  (2.0d0*pi) 

i 

!  Calculate  synthesized  tine  history  by  inverse  FFT 

i 

CALL  DFFTCB  (ran,  freg,  wregm) 
wregr=0 . 0 
do  i=2,mm 

ii=mm-i+l 

wregr  (i)  =dble  (wregm  (i) ) 
wregi  (i)  =dimag  (wregm  (i) ) 
enddo 

i 

!  Output  synthesized  time  histories 

i 

do  i=2,mm/2 

write  (6, 200)  tt(i),  wregr  (i) ,  wregi  (i) ,  w(i) 

200  format (5 (lx,el5.7) ) 

enddo 

i 

!  Sort  and  Output  all  modal  parameters 

! 

call  sortmodal  (nq,  wr,  xr,  cmodal,  arn) 
do  i=l,nq 

write  (4, 100)  i,  wr  (i)  /  (2 . 0d0*pi) ,  xr(i),  cmodal  (i) ,  arn(i) 
100  format (li5,  5 (lx,  el4.6)) 

enddo 

close (unit=3) 
close (unit=4) 
close (unit=6) 

8888  continue 

877  format (5 (lx, el5. 7)) 

400  format (5 (lx,el5.7) ) 
close (unit=l) 
close (unit=2) 
close (unit=5) 

stop 

end 


Subroutine  to  calculate  damping  curve  and  standard  deviation 


Subroutine  meansd(nto,  totalfreq,  totaldamp,  rayleigh,  raywgt,  & 
ffreq,  fdampl,  fdamp2,  sdray,  sdwgt) 
implicit  double  precision  (a-h,o-z) 

dimension  totalfreq (2000) ,  totaldamp (2000) ,  ffreq (250),  fdampl (250) ,  & 
fdamp2  (250) 

dimension  rayleigh ( 2 ) ,  raywgt ( 2 ) 
pi=4 . 0d0*datan (1 . OdO) 
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sdray=0 .  OdO 
sdwgt=0 .  OdO 
do  i=l,nto 

w=2 . OdO*pi*totalfreq (i) 

x=rayleigh (1) /2 . OdO/w+rayleigh (2) /2 . OdO*w 
sdray=sdray+ (totaldamp (i) -x) **2 . OdO 
x=raywgt ( 1 ) / 2 . OdO/w+raywgt (2 ) / 2 . OdO*w 
sdwgt=sdwgt+ (totaldamp (i) -x) **2 . OdO 
enddo 

sdray=sdray/nto 
sdwgt=sdwgt/nto 
do  i=l,250 
f=dfloat (i) 
ffreq(i)  =f 
w=2 . OdO*pi*f 

fdampl (i) =rayleigh (1) /2 . OdO/w+rayleigh (2) /2 . OdO*w 
fdamp2 (i) =raywgt (1) /2 . OdO/w+raywgt (2) /2 . OdO*w 
enddo 
return 
end 


!  Sort  total 


subroutine  sorttotal (itt,  ffreq,  ddamp,  ar) 
implicit  double  precision  (a-h,o-z) 
dimension  ffreq(2000),  ddamp(2000)  ,  ar(2000) 
ntt=itt 
do  i=l,ntt-l 
do  j=i+l,  ntt 

if  (ffreq (i)  .gt. ffreq (j) )  then 
ddum=ffreq (i) 

ffreq (i) =ffreq ( j ) 
ffreq ( j ) =ddum 
ddum=ddamp  (i) 
ddamp  (i)  =ddamp  ( j ) 
ddamp  ( j )  =ddum 
ddum=ar (i) 
ar  (i)  =ar  (j ) 
ar  (j)=ddum 

endif 

enddo 

enddo 

return 

end 

i 

i - 

subroutine  sortmodal  (nq,  wr,  xr,  cmodal,  arn) 


implicit  double  precision  (a-h,o-z) 
dimension  wr (2000) ,  xr(2000)  ,  arn (2000) 
double  complex  cmodal (2000) ,  cdum 
do  i=l,nq-l 
do  j=i+l,  nq 

if  (wr  (i)  .gt  .wr  (j ) )  then 
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ddum=wr  (i) 

wr  (i)  =wr  (j ) 
wr  (j )  =ddum 
ddum=xr (i) 
xr  (i)  =xr  (j ) 
xr ( j ) =ddum 
ddum=arn (i) 
arn  (i)  =arn  (j ) 
arn ( j ) =ddum 
cdum=cmodal ( i ) 
cmodal  (i)  =cmodal  ( j ) 
cmodal  ( j )  =cdum 

endif 

enddo 

enddo 

return 

end 

i 

i - 

subroutine  FindAr  (nq,  ran,  dt ,  w,  sr ,  cmodal ) 


!  to  find  modal  constants  with  frequency  domain  least  square  method 
implicit  double  precision  (a-h, o-z) 
dimension  w(2000) 

double  complex  sr (201) , cmodal (201) ,t (2000) ,ss (2000,200) ,  & 

st  (200, 2000) ,  stf  (200) ,  sts  (200,200) ,  ff  (2000) , wi 
pi=4 . 0d0*datan  (1 .  OdO) 
do  i=l,mm 

t  (i)  =dcmplx  (w  (i) ,  0 .  OdO) 
enddo 

CALL  DFFTCF  (ran,  t,  ff) 
m=mm/2 

df=l .  0d0/  (dt*dfloat  (ran-1) ) 
dw=2 .  0d0*pi*df 
ff  (1  :mm)  -ff  (1  :ran)  /dfloat  (mm) 
do  felL, mil 

w=dfloat (i— 1) *dw 
wi=dcmplx  ( 0 .  OdO ,  w) 
do  j=l,nq 

ss  (i,  j )  =wi/  (wi-sr  ( j ) )  *df 
enddo 

enddo 

do  i=mm,mi2,-l 

w=  dfloat  (mm-iil)  *dw 
wi=dcnplx  (0 .  OdO,  -w) 
do  j=l,nq 

ss  (i,  j )  =wi/  (wi-sr  ( j ) )  *df 
enddo 

enddo 

st  (1  :nq,  1  :mm)  =transpose  (ss  (1  :mm,  1  :nq) ) 

stf  (1  :nq)  ^matmul  (st  (1  :nq,  1  :mm) ,  ff  (1  :rrm) ) 

sts  (1  :nq,  1  :nq)  ^matmul  (st  (1  :nq,  1  :mm) ,  ss  (1  :mm,  l:nq) ) 

CALL  DLSACG  (nq,  sts,  100,  stf,  1,  cmodal) 
return 
end 

i 

i - 
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SUBROUTINE  BAUPOL  (COEERE,  COEFIM, N,  COMPL,  ROOTRE,  ROOTIM,  NUMIT) 

I  Q*  *************************************************************************  •*••*••*• 

!C  * 

!C  Without  knowing  any  approximations  of  the  roots,  this  * 

!C  SUBROUTINE  finds  N  approximate  values  Z(I),  1=1,  ...,  N  for  * 

!C  the  N  zeros  of  a  polynomial  PN  of  degree  N  with  real  or  * 

!C  complex  coefficients.  * 

!C  The  polynomial  is  described  as  follows:  * 

!C  * 

!C  PN (Z) =COEF (0) +COEF (1) *Z+COEF (2) *Z**2+. . . +COEF (N) *Z**N,  * 

!C  * 

!C  with  COEF(I)  =  (COEERE  ( I ) ,  COEFIM  ( I ) )  for  1=0,  ...,  N  (complex  * 

!C  coefficients) .  * 

!C  * 

!C  INPUT  PARAMETERS:  * 

!C  ===========  * 

!C  COEERE  :  (N+l) -vector  COEERE(0:N)  containing  the  real  * 

!C  part  of  each  coefficient  of  the  polynomial  PN  in  * 

!C  DOUBLE  PRECISION.  * 

!C  COEFIM  :  (N+l) -vector  COEFIM  (0:N)  containing  the  imaginary  * 

!C  part  of  each  coefficient  of  the  polynomial  PN  in  * 

!C  DOUBLE  PRECISION.  * 

!C  N  :  degree  of  the  polynomial  PN.  * 

!C  COMPL  :  boolean  variable  :  * 

!C  COMPI^.TRUE.  ,  if  the  coefficients  are  COMPLEX.  * 

!C  COMPI^ . FALSE .  ,  if  the  coefficients  are  REAL.  * 

!C  * 

!C  OUTPUT  PARAMETERS:  * 

!C  ======  * 

!C  ROOTRE  :  N-vector  ROOTRE(l:N)  containing  the  approximate  * 

!C  real  parts  of  the  computed  zeros  of  PN  in  * 

!C  DOUBLE  PRECISION.  * 

!C  ROOTIM  :  N-vector  ROOTIM(l:N)  containing  the  approximate  * 

!C  imaginary  parts  of  the  computed  zeros  of  PN  in  * 

!C  DOUBLE  PRECISION.  * 

!C  NUMIT  :  maximum  number  of  iteration  steps.  * 

!C  * 

IMPLICIT  DOUBLE  PRECISION  (a-h,o-z) 

DIMENSION  COEERE  (0:201) ,  COEFIM  (0:201)  ,E  (201) ,  & 

ROOTRE (200) ,ROOTIM(200) ,A(500) ,B (500) , C (500) 
DOUBLE  PRECISION  INFINY,  MAX,  MIN 
LOGICAL  COMPL 

i 

!  Initializing  the  iteration  step  counter  NUMIT  and  the  error 
!  bound  GAMMA. 

i 

NUMIT=0 
GAMMA=5 . OD-18 
IF (N  .EQ.  1)  THEN 

i 

!  If  the  degree  of  PN  is  N=l,  then  the  zero  of  the  polynomial  PN  is 
!  Z=-COEF ( 0 ) / COEF ( 1 ) ,  where  COEF (I) = (COEERE (I) , COEFIM (I) )  for  1=0,1. 

i 

CALL  CDIV  (-COEERE  (0)  ,  -COEFIM  (0)  ,  COEERE  (1)  ,  COEFIM  (1)  ,  & 
ROOTRE  (1)  ,  ROOTIM  (1)  ) 

RETURN 
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ELSE 

N1=N+1 

| 

!  Scaling  via  SCALFC . 

j 

DO  10  1=1, N1 

E  (I)  =ABSCOM  (COEFRE  (Nl-I) ,  COEFIM  (Nl-I)  ) 

10  CONTINUE 

CALL  MCONST (EMACHP, INFINY, SMALNO, BASE) 

BND=SCALFC (N1,E,  EMACHP, INFINY, SMALNO, BASE) 

IF (BND  .EQ.  1.0D0)  THEN 

| 

!  Normalizing,  in  case  scaling  by  SCALFC  did  not  normalize  the  coefficients. 

I 

MAX=0 . 0D0 
MIN=1 . 0DH-300 
DO  20  I=N, 0,  -1 

X=ABSCOM  (COEFRE  (I)  ,  COEFIM  (I)  ) 

IF (X  .GT.  MAX)  MAX=X 

IF (X  .LT.  MIN  .AND.  X  .NE.  0.0D0)  MIN=X 
20  CONTINUE 

BND=DSQRT  (MAX*MIN) 

BND=1 . 0D0/BND 
END  IF 
DO  30  K=1,N1 

B  (2*K-1)  COEFRE  (Nl-K)  *BND 
B  (2*K)  =0 . 0D0 

IF  (CCMPL)  B  (2*K)  COEFIM  (Nl-K)  *BND 
30  CONTINUE 
X0=0 . 0D0 
Y0=0 . 0D0 
DO  40  1=1, N 
L=2*  (N+2-I) 

DO  50  K=1,L 
A  (K)  =B  (K) 

50  CONTINUE 

| 

!  Calculating  the  I-th  zero  of  PN. 

I 

CALL  BAUZRO (X0, Y0, N+l-I, GAMMA, XNEW, YNEW,NUMIT, A, B, C) 

ROOTRE  (I)  =XNEW 
ROOTIM  ( I )  =YNEW 
X0=XNEW 
Y0=-YNEW 
40  CONTINUE 
END  IF 
RETURN 
END 


SUBROUTINE  BAUZRO (X0 , Y0, N, GAMMA, XNEW, YNEW, NUMIT, A, B, C) 

* 

This  SUBROUTINE  calculates  a  zero  of  a  polynomial  PN  with  * 
complex  coefficients .  * 

We  solve  the  equation  PN (Z) /PN' (Z) =0  * 
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via  Newton's  method  with  spiralization  and  extrapolation  to  * 
improve  convergence .  * 
The  initial  approximation  (X0+I*Y0)  can  be  chosen  arbitrarily.* 

* 


INPUT  PARAMETERS: 

* 

* 

XO 

real  part  of  the  initial  approximation. 

* 

Y0 

imaginary  part  of  the  initial  approximation. 

* 

N 

degree  of  the  polynomial  PN. 

* 

GAMMA 

error  bound. 

* 

* 

OUTPUT  PARAMETERS: 

* 

XNEW 

real  part  of  the  computed  zero  of  PN. 

* 

YNEW 

imaginary  part  of  the  zero  of  PN. 

* 

NUMIT 

maximum  number  of  iteration  steps. 

* 

* 


IMPLICIT  DOUBLE  PRECISION  (a-h,o-z) 
DIMENSION  A(500) ,B (500) ,C (500) 
LOGICAL  ENDIT 

IF (N  .EQ.  2)  THEN 


Solving  the  remaining  2nd  degree  polynomial  exactly. 

CALL  CDIV (A (3)  ,A(4)  ,A(1)  , A(2) , P1RE, P1IM) 

CALL  CDIV  (A  (5)  ,A(6)  ,A(1)  ,A(2)  ,Q1RE,Q1IM) 
P12RE=-PlRE/2 . 0D0 
P12IM=-PlIM/2 . 0D0 
RA1RE=P12RE*P12EE-P12IM*P12IM 
RA1IM=2 . 0D0*P12RE*P12IM 
RARE=RA1RE-Q1RE 
RAIM=RA1 IM-Q1 IM 
IF (RAIM  .EQ.  0.0D0)  THEN 
IF (RARE  .LT.  0.0D0)  THEN 


Purely  imaginary  root. 

RTIM=DSQRT (-RARE) 
RTRE=0 . 0D0 
XNEW=P12RE 
YNEW=P12  IM+RTIM 
RETURN 
ELSE 


Real  root. 


RTRE=DSQRT (RARE) 
RTIM=0 . 0D0 
XNEW=P12RE+RTRE 
YNEW=P12IM 
RETURN 
END  IF 
ELSE 
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Conplex  root. 


RABE=ABSCOM (RARE, RAIM) 

IF (RARE  .GT.  0.0D0)  THEN 

RTRE=DSQRT (0 . 5D0* (RABE+RARE) ) 
IF (RAIM  .LT.  0.0D0)  RTRE=-RTRE 
RTIM=0 . 5D0  *RAIM/RTRE 
ELSE 

RTIM=DSQRT (0 . 5D0* (RABE-RARE) ) 
RTRE=0 . 5D0*RAIM/RTIM 
END  IF 

XNEW=P1 2RE+RTRE 
YNEW=P12IM+RTIM 
RETURN 
END  IF 

ELSE  IF (N  .EQ.  1)  THEN 


Polynomial  of  1st  degree. 

XNEW=P1 2RE-RTRE 
YNEW=P1 2 IM-RT IM 
RETURN 
ELSE 
1=0 

ENDIT=. FALSE. 

RHODSQRT  (GAMMA) 

BETA=10 . 0D0*GAMMA 
QR=0 . 1D0 
QI=0 . 9D0 
XNEW=X0 
YNEW=Y0 

CALL  COMHOR  (XNEW,  YNEW,  N,  GAMMA,  UNEW,  VNEW,  UDNEW,  VDNEW,  & 
UDDNEW,  VDDNEW,  BD,  BDD,  A,  B,  C) 

NUMIT=NUMIT+1 
PBNEW=ABSCCM  (UNEW,  VNEW) 

IF (PBNEW  .LE.  BD)  THEN 
RETURN 
ELSE 

PBOLD=2 . 0D0*PBNEW 

DZMIN=BETA*  (RHO+ABSCOM  (XNEW,  YNEW)  ) 

1 0  PSBNEW=ABSCQM  (UDNEW,  VDNEW) 

Spiralization . 

IF (PBNEW  .LT.  PBOLD)  THEN 

DZMAX=1 . 0D0+ABSCCM  (XNEW,  YNEW) 

NUMIT=NUMIT+1 

H1=UDNEW*UDNEW-VDNEW*VDNEW-UNEW*UDDNEW+VNEW*VDDNEW 
H2=2 . ODO*UDNEW*VDNEW-UNEW*VDDNEW-VNEW*UDDNEW 
IF (PSBNEW  .GT.  10.0D0*BDD  .AND.  & 

ABSCOM (HI , H2 )  .GT.  100 . 0D0*BDD*BDD)  THEN 


Applying  Newton's  method. 


1=1+1 

IF (I  .GT.  2)  1=2 
U=UNEW*UDNEW-VNEW*VDNEW 
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V=UNEW^VDNEW+VNEW*UDNEW 
CALL  CDIV  (-U,  -V,  HI ,  H2 ,  DX,  DY) 

IF  (ABSCOM  (DX,  DY)  .GT.  DZMAX)  THEN 
H=DZMAX/ABSCOM  (DX,  DY) 

DX=DX*H 
DY=DY*H 
1=0 
END  IF 

IF  (I  .EQ.  2  .AND.  ABSCOM  (DX,  DY)  .LT.  DZMIN/RHO  .AND.  & 
ABSCOM (DX,DY)  .GT.  0.0D0)  THEN 

Extrapolation . 

1=0 

CALL  CDIV  (XNEW-XOLD,  YNEW-YOLD,  DX,  DY,  H3 ,  H4 ) 

H3=l . 0D0+H3 

Hl=H3*H3-H4*H4 

H2=2.0D0*H3*H4 

CALL  CDIV (DX, DY, Hi , H2 , H3 , H4 ) 

IF  (ABSCOM (H3,H4)  .LT.  50 . 0D0*DZMIN)  THEN 
DX=DX+H3 
DY=DY+H4 
END  IF 
END  IF 
XOLD=XNEW 
YOLD=YNEW 
PBOLD=PBNEW 
ELSE 


In  a  neighborhood  of  a  saddle  point. 


1=0 

H=DZMAX/ PBNEW 

DX=H*UNEW 

DY=H*VNEW 

XOLD=XNEW 

YOLD=YNEW 

PBOLD=PBNEW 

20  CALL  COMHOR  (XNEW+DX,  YNEW+DY,  N,  GAMMA,  U,  V,  H,  Hi ,  & 

H2,H3,H4,H5,A,B,C) 

IF  (DABS  (ABSCOM  (U,  V) /PBNEW-1 . 0D0)  .LE.  RHO)  THEN 
DX=2 . 0D0*DX 
DY=2.0D0*DY 
GOTO  20 
END  IF 
END  IF 
ELSE 
1=0 

NUMI  T=NUMI  T+ 1 
H=QR*DX-QI *DY 
DY=QR*  DY+QI * DX 
DX=H 
END  IF 

IF (ENDIT)  THEN 

IF (ABSCOM (DX,DY)  .LT.  0.1D0*BDZE)  THEN 
XNEW=XNEW+DX 
YNEW=YNEW+DY 
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END  IF 

CALL  COMHOR  (XNEW,  YNEW,  N,  GAMMA,  UNEW,  VNEW,  UDNEW,  & 
VDNEW,  UDDNEW,  VDDNEW,  BD,  BDD,  A,B,C) 

RETURN 

ELSE 

XNEW=XOLD+DX 

YNEW=YOLD+DY 

DZMIN=BETA*  (RHO+ABSCOM  (XNEW,  YNEW)  ) 

CALL  COMHOR  (XNEW,  YNEW,  N,  GAMMA,  UNEW,  VNEW,  UDNEW,  & 
VDNEW,  UDDNEW,  VDDNEW,  BD,  BDD,  A,  B,  C) 
PBNEW=ABSCOM  (UNEW,  VNEW) 

IF (PBNEW  .EQ.  0.0D0)  THEN 
RETURN 

ELSE  IF  (ABSCOM (DX,  DY)  .GT.  DZMIN  .AND.  & 

PBNEW  .GT.  BD)  THEN 
GOTO  10 
ELSE 

ENDIT= . TRUE . 

BDZE=ABSCOM (DX, DY) 

GOTO  10 
END  IF 
END  IF 
END  IF 
END  IF 
END 


SUBROUTINE  COMHOR  (X,  Y,  N,  GAMMA,  U, V,  UD,  VD,  UDD, VDD,  BDP,  BDPD, A,  B,  C) 


i?ititiritit-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

~k 


!  This  SUBROUTINE  calculates  the  complex  functional  value  * 
!  PN  (X+I*Y)  =  (U+I*V) ,  the  complex  valued  derivatives  * 
!  PN'  (X+I*Y)  =  (UD+I*VD)  and  PN'  '  (X+I*Y)  =  (UDD+I*VDD)  of  a  polymial* 
!  PN(Z)  of  degree  N  (N>0)  with  complex  coefficients  by  using  * 
!  the  Horner-scheme .  * 
!  Additionally  bounds  BDP  and  BDPD  are  computed  for  the  * 
!  rounding  errors  in  computing  DABS(PN(Z))  and  DABS(PN' (Z) ) .  * 
!  The  complex  coefficients  of  PN  are  stored  in  a  2-dimensional  * 
!  array  A (1:2* (N+l) ) ,  arranged  in  descending  order  of  the  * 
!  powers  (they  will  remain  unchanged  by  this  subroutine) .  * 
!  The  complex  coefficients  of  the  polynomial  Q(Z)  of  degree  N-l  * 
!  are  stored  in  the  array  B (1:2* (N+l) ) ,  arranged  in  descending  * 
!  order.  Here  Q(Z)  is  defined  by  PN (Z)  =Q (Z)  *  (Z-ZO)  +PN (Z0)  .  * 
!  The  array  C(1:2*N)  is  used  as  an  auxiliary  array.  * 

i  * 
i  * 


!  INPUT  PARAMETERS: 


!  X  :  real  part  of  the  value  for  which  the  functional 

!  value  and  its  1st  and  2nd  derivatives  are  to  be 

!  computed  for  the  polynomial  PN. 

!  Y  :  imaginary  part  of  the  value  for  which  the 

!  functional  value  and  its  1st  and  2nd  derivatives 

!  are  to  be  computed  for  the  polynomial  PN. 

!  GAMMA  :  error  bound. 


* 

* 

* 

* 

* 

* 

* 

* 

* 
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OUTPUT  PARAMETERS: 


* 


* 

* 


u 

:  real  part  of  PN (X+I*Y) . 

* 

V 

:  imaginary  part  of  PN (X+I*Y) . 

* 

* 

UD 

:  real  part  of  PN' (X+I*Y) . 

* 

VD 

:  imaginary  part  of  PN' (X+I*Y) . 

* 

* 

UDD 

:  real  part  of  PNf 1 (X+I*Y) . 

* 

VDD 

:  imaginary  part  of  PN'  ' (X+I*Y) . 

* 

* 

BDP 

:  bound  for  the  rounding  error  of 

DABS  (PN  (Z)  )  . 

* 

BDPD 

:  bound  for  the  rounding  error  of 

DABS (PN'  (Z)) . 

* 

* 


IMPLICIT  DOUBLE  PRECISION  (a-h,o-z) 
DIMENSION  A(500) ,B (500)  ,C (500) 

C  (1)  =A  (1) 

B  (1)  =A  (1) 

C  (2)  =A(2) 

B  (2)  =A(2) 

BDPD=ABSCOM(A(l) ,A(2) ) 

BDP=BDPD 

MS=N-1 

M=N 

J=N 

NM2P1=N*2+1 
DO  10  K=3,NM2P1,2 
J=J-1 

H1=X*B (K-2) -Y*B (K-l) 

H2=Y*B (K-2) +X*B (K-l) 

B  (K)  =A  (K)  +H1 
B  (K+l )  =A  (K+l )  +H2 
H3=ABSCCM  (A  (K) ,  A  (K+l ) ) 

H4=ABSCCM (HI , H2 ) 

H=H3 

IF(H3  .LT.  H4)  H=H4 
IF (H  .GT.  BDP)  THEN 
BDP=H 
MNJ 
END  IF 

IF (K  .EQ.  NM2P1)  THEN 
GOTO  20 
ELSE 

H1=X*C (K-2) -Y*C (K-l) 

H2=Y*C (K-2) +X*C (K-l) 

C  (K)  =B  (K)  +H1 
C  (K+l)  =B  (K+l)  +H2 
H3=ABSCCM (B  (K) ,  B  (K+l)  ) 
H4=ABSCCM  (HI , H2 ) 

H=H3 

IF(H3  .LT.  H4)  H=H4 
IF (BDPD  .LT.  H)  THEN 
BDPD=H 
MS=J-1 
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END  IF 
END  IF 
10  CONTINUE 
20  CONTINUE 
U=B (2*N+1) 

V=B (2*N+2) 

UD=C(2*N-1) 

VD=C  (2*N) 

H=ABSCOM  (X,  Y) 

IF (H  .NE.  0.0D0)  THEN 

BDP=BDP* FLOAT (M+l) *H**M 
BDPI>BDPD*FLOAT (MS+1) *H**MS 
ELSE 

BDP=ABSCCM(U,V) 

BDPD=ABSCOM (UD,  VD) 

END  IF 

BDP=BDP*GAMYIA 
BDPD=BDPD*GAMMA 
IF (N  .GT.  1)  THEN 
H1=C  (1) 

H2=C  (2) 

NM2M3=N*2-3 
DO  30  K=3,NM2M3, 2 
H=C (K) +X*H1-Y*H2 
H2=C (K+l) +Y*H1+X*H2 
H1=H 

30  CONTINUE 

UDD=2.0D0*H1 

VDD=2.0D0*H2 

RETURN 

ELSE 

UDD=0 . 0D0 
VDD=0 . 0D0 
RETURN 
END  IF 
END 

| 

i 


DOUBLE  PRECISION  FUNCTION  ABSCCM(X,Y) 

* 

This  FUNCTION- subroutine  calculates  the  absolute  value  of  a  * 
complex  number  (X+I*Y)  .  * 

* 


INPUT  PARAMETERS:  * 

========  * 

X  :  real  part  of  the  complex  number.  * 

Y  :  imaginary  part  of  the  complex  number.  * 

* 


OUTPUT  PARAMETER:  * 

========  * 

ABSCOM  :  absolute  value  of  the  complex  number.  * 

* 
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IMPLICIT  DOUBLE  PRECISION  (a-h,o-z) 

IF (X  .NE.  0.0D0  .OR.  Y  .NE.  0.0D0)  THEN 
IF (DABS (X)  .GE.  DABS (Y) )  THEN 

ABSCOM=DABS (X) *DSQRT (Y/X*Y/X+1 . 0D0) 
RETURN 
ELSE 

ABSCOM=DABS  (Y)  *DSQRT  (X/Y*X/Y+1 . 0D0) 
RETURN 
END  IF 
ELSE 

ABSCOM=0 . 0D0 
RETURN 
END  IF 
END 


SUBROUTINE  MCONST  (FMACHP,  INFINY,  SMAINO,  BASE) 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

k 

This  subroutine  sets  up  some  constants  that  are  machine  * 

dependent .  * 

k 


k 


OUTPUT  PARAMETERS:  * 

-  k 


FMACHP 

INFINY 

SMALNO 

BASE 


machine  constant  for  DOUBLE  PRECISION.  * 
largest  representable  floating-point  number.  * 
smallest  representable  floating-point  number.  * 
base  of  the  floating-point  number  system  used  to  * 
represent  machine  numbers.  * 

k 


k 

Description  of  the  auxiliary  variables:  * 

-  k 


I 

M 

N 


:  number  of  digits  of  the  floating-point  mantissa  * 
of  DOUBLE  PRECISION  numbers.  * 

:  largest  allowed  exponent.  * 

:  smallest  allowed  exponent.  * 

k 


IMPLICIT  DOUBLE  PRECISION  (a-h,o-z) 

DOUBLE  PRECISION  INFINY 
BASE=8 . 0D0 
1=29 
M=322 
N=-2 93 
FMACHP=1 . 0D0 
10  FMACHP=0 . 5D0*EMACHP 

IF (1 . 0D0  .LT.  1 . 0D0+FMACHP)  GOTO  10 
FMACHP=2 . 0D0*EMACHP 

INFINY=BASE*  (1 .  0D0-BASE**  (-1)  )  *BASE**  (M-l) 

SMALNO=  (BASE**  (N+3)  )  /BASE**3 

RETURN 

END 
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I 


DOUBLE  PRECISION  FUNCTION  SCALFC (NN, PT, FMACHP, INFINY, SMALNO,  & 

BASE) 

i?itititiriti?itit-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

~k 

This  FUNCTION- subroutine  calculates  a  scaling  factor  which  * 
is  used  to  scale  the  polynomial  coefficients.  * 

* 

* 

INPUT  PARAMETERS:  * 


NN 

PT 

FMACHP 

INFINY 

SMALNO 

BASE 


1  +  the  degree  of  the  polynomial, 
nn-vector  PT(1:NN)  containing  the  absolute 
values  of  the  polynomial's  coefficients, 
machine  constant  for  DOUBLE  PRECISION, 
largest  representable  floating-point  number, 
smallest  representable  floating-point  number, 
base  for  the  floating-point  number  system  used  by 
the  machine. 


OUTPUT  PARAMETER: 


SCALFC  :  scaling  factor. 


IMPLICIT  DOUBLE  PRECISION  (a-h,o-z) 
DIMENSION  PT  (NN) 

DOUBLE  PRECISION  MAX,  MIN  ,  INFINY 
HI=DSQRT  (INFINY) 

LO=SMALNO/ EMACHP 
MAX=0 . 0D0 
MIN=INFINY 
DO  10  1=1, NN 
X=PT(I) 

IF (X  .GT.  MAX)  MAX=X 

IF (X  .NE.  0.0D0  .AND.  X  .LT.  MIN)  MIN=X 
10  CONTINUE 
SCALFC=1 . 0D0 

IF  (MIN  .GE.  LO  .AND.  MAX  .LE.  HI)  THEN 
RETURN 
ELSE 

X=LO/MIN 

IF (X  .GT.  1.0D0)  THEN 
SC=X 

IF ( INFINY/ SC  .GT.  MAX)  SC=1.0D0 
ELSE 

SC=1 . 0D0/  (DSQRT  (MAX)  *DSQRT  (MIN)  ) 

END  IF 

L=DLOG (SC) /DLOG (BASE) +0 . 5D0 
SCALFC=BASE**L 
END  IF 
RETURN 
END 
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SUBROUTINE  CDIV(A,B,C,D,X,Y) 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* 

This  SUBROUTINE  performs  a  complex  division  * 


(X+I*Y)  :=  (A+I*B)  /  (C+I*D)  . 

* 

k 

INPUT 

PARAMETERS : 

k 

k 

k 

A 

:  real  part  of  the  numerator. 

k 

B 

:  imaginary  part  of  the  numerator. 

k 

k 

C 

:  real  part  of  the  denominator. 

k 

D 

:  imaginary  part  of  the  denominator. 

k 

k 

OUTPUT  PARAMETERS: 

k 

k 

k 

X 

:  real  part  of  the  quotient. 

k 

Y 

:  imaginary  part  of  the  quotient. 

k 

* 

* 


NOTE:  If  the  denominator's  real  and  imaginary  parts  are  both  * 
equal  to  zero,  the  program  is  aborted  with  a  detailed  * 

error  message .  * 

* 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


IMPLICIT  DOUBLE  PRECISION  (a-h,o-z) 

IF (C  .NE.  0.0D0  .OR.  D  .NE.  0.0D0)  THEN 
IF (A  .NE.  0.0D0  .OR.  B  .NE.  0.0D0)  THEN 
IF  (DABS  (A)  .GT.  DABS  (B)  )  THEN 
U=A 

AM=1 . 0D0 
AN=B/A 
ELSE 
U=B 

AM=A/B 
AN=1 . 0D0 
END  IF 

IF  (DABS  (C)  .GT.  DABS  (D)  )  THEN 
V=C 

P=1 . 0D0 
Q=D/C 
ELSE 
V=D 
P=C/D 
0=1 . 0D0 
END  IF 
F=U/V 
V=P*P+Q*Q 
U= (AM*P+AN*Q) /V 
X=U*F 

U= (-AM*Q+AN*P) /V 
Y=U*F 
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RETURN 

ELSE 

X=0 . ODO 
Y=0 . ODO 
RETURN 
END  IF 
ELSE 

WRITE  (*,  *)  'DIVISION  BY  ZERO  IN  SUBROUTINE  CDIV' 
STOP 
END  IF 
END 
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